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We investigate non-equilibrium steady states of driven-dissipative ideal quantum gases of both 
bosons and fermions. We focus on systems of sharp particle number that are driven out of equilibrium 
either by the coupling to several heat baths of different temperature or by time-periodic driving in 
combination with the coupling to a heat bath. Within the framework of (Floquet-)Born-Markov 
theory, several analytical and numerical methods are described in detail. This includes a mean- 
held theory in terms of occupation numbers, an augmented mean-held theory taking into account 
also non-trivial two-particle correlations, and quantum-jump-type Monte-Carlo simulations. For the 
case of the ideal Fermi gas, these methods are applied to simple lattice models and the possibility of 
achieving exotic states via bath engineering is pointed out. The largest part of this work is devoted 
to bosonic quantum gases and the phenomenon of Bose selection, a non-equilibrium generalization of 
Bose condensation, where multiple single-particle states are selected to acquire a large occupation 
[Phys. Rev. Lett. Ill, 240405 (2013)]. In this context, among others, we provide a theory for 
transitions where the set of selected states changes, describe an efficient algorithm for hnding the 
set of selected states, investigate beyond-mean-held effects, and identify the dominant mechanisms 
for heat transport in the Bose selected state. 

PACS numbers: 05.30.Jp, 05.70.Ln, 67.10.Ba, 67.85.Jk 


I. INTRODUCTION 

There is a huge current interest in non-equilibrium phe¬ 
nomena of many-body systems beyond the hydrodynamic 
description of systems retaining approximate local equi¬ 
librium. Recent work concerns several paradigmatic sce¬ 
narios, like the dynamics away from equilibrium in re¬ 
sponse to a slow or an abrupt parameter variation 
the possible relaxation towards equilibrium HE! versus 
many-body localization ma, and the control of many- 
body physics by means of strong periodic forcing [MI]. 
Also the possibility to achieve transient light-induced su¬ 
perconductivity above the equilibrium critical tempera¬ 
ture attracted enormous interest m- 

Another fundamental scenario of many-body dynamics 
are driven-dissipative quantum systems and their non¬ 
equilibrium steady states [TBfl22j . These include, for ex¬ 
ample, time-periodically driven open many-body systems 
[23H27] and photonic many-body systems [28H33] . In 
contrast to equilibrium states, which depend on a few 
thermodynamic parameters like temperature and chemi¬ 
cal potential only, such non-equilibrium steady states de¬ 
pend on the very details of the environment. On the one 
hand, this makes their theoretical treatment challenging. 
On the other hand, it offers also interesting opportunities 
to engineer the state and the properties of a many-body 
system beyond the constraints of thermal equilibrium in 
a robust and controlled fashion. 

In this context, it was recently pointed out that already 
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FIG. 1. (color online) Two paradigmatic examples of driven- 
dissipative ideal quantum gases possessing non-equilibrium 
steady states, (a) Periodically driven system weakly coupled 
to a heat bath, (b) Autonomous system weakly coupled to 
two heat baths of different temperature. 


an ideal Bose gas of N particles can exhibit intriguing be¬ 
havior, when it is driven into a steady state far from equi¬ 
librium, e.g., by coupling it to two heat baths of different 
temperature or by time-periodic driving in the presence 
of a heat bath (see Fig. [^. In the quantum degener¬ 
ate regime of large densities, the Bose gas undergoes a 
generalized form of Bose condensation, where multiple 
single-particle states can be selected to acquire large oc¬ 
cupations [53]. Namely, the single-particle states unam¬ 
biguously separate into two groups: one that is called 
Bose selected, whose occupations increase linearly when 
the total particle number is increased at fixed system 
size, and another one whose occupations saturate. This 
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phenomenon is a consequence of the bosonic quantum 
statistics. It includes standard Bose condensation into a 
single quantum state, fragmented Bose condensation into 
a small number of single-particle states each acquiring a 
macroscopic occupation, and the case where a fraction of 
all single-particle states acquires large, but individually 
non-extensive occupations. The properties of the system, 
like its coherence or its heat conductivity, sensitively de¬ 
pend on which of these scenarios occurs. 

The physics of driven-dissipative ideal Bose gases is 
intimately related also to collective effects in classical 
systems and processes, where bunching phenomena have 
been identified as analog of Bose condensation. This in¬ 
cludes the dynamics of networks and economic models 
[331 |5S], classical transport and traffic chemi¬ 

cal reactions [45], as well as population dynamics and 
evolutionary game theory |46j . These connections have 
recently been discussed nicely by Knebel et al. m- 


In this paper, we investigate non-equilibrium steady 
states of driven-dissipative ideal quantum gases of both 
bosons and fermions. We focus on systems of sharp parti¬ 
cle number that exchange energy with the environment. 
These quantum gases are driven out of equilibrium ei¬ 
ther by the coupling to several heat baths of different 
temperature or by time-periodic driving in combination 
with the coupling to a heat bath (see Fig. [^. We treat 
the problem using (Floquet-)Born-Markov theory [48l - 
I52j . which is valid in the limit of weak system-bath cou¬ 
pling. In section jn] this theoretical framework is reviewed 
and applied to the problem of the ideal quantum gas. 
Morever, several model systems are introduced. In or¬ 
der to treat the resulting many-body master equation, 
we then describe analytical and numerical methods for 
computing the steady state (Section III). This includes 
a standard mean-field description in terms of single¬ 
particle occupation numbers. We, moreover, derive an 
augmented mean-field theory taking into account also 
non-trivial two-particle correlations, and explain how to 
apply quantum-jump-type Monte-Carlo simulations to 
the problem. These methods are then applied to both 
the ideal Bose gas (Section IV) and the ideal Fermi gas 
(Section [vj). 

Our treatment of the fermionic case in Section El is 
rather brief and demonstrates the application of our the¬ 
ory to simple lattice models and the possibility to achieve 
exotic states via bath engineering. These results can 
be relevant, e.g., for the problem of realizing Floquet 
topological insulators with periodically forced electronic 
systems (graphene |7] or semiconductor heterostructures 

m)- 


The largest part of this paper is devoted to bosonic 
quantum gases and the phenomenon of Bose selection 
discussed in Section jlV] Here we first review equilibrium 
Bose condensation (Sec. IV AI and Bose selection in non¬ 
equilibrium steady states (Secs. IV B to IV Ej gi ve a de¬ 
tailed discussion of the results of Reference 323] 1. After 
that, we derive a theory for transitions where the set of 
selected states changes (Sec. IV F), present an efficient al¬ 


gorithm for finding the set of selected states (Sec. IV G I, 
discuss the possibility of approaching a preasymptotic 
state at intermediate densities before the tr ue asy mp- 
totic state is reached at large densities (Sec. IVHi, in¬ 
vestigate the properties of sys tems described by non-fully 
connected rate matrices (Sec. IV11, study the role o f fluc¬ 
tuations and beyond mean-field effects (Sec. IV J), and 
identify the dominant mechanisms for heat transport in 
the Bose selected state Sec. nvij 


II. GENERAL FRAMEWORK AND MODELS 

In this section we set up the master equations for 
an ideal quantum gas of N indistinguishable, nonin¬ 
teracting particles, weakly coupled to one or several 
heat baths. We cover both the case of an autonomous 
system with time-independent Hamiltonian H and the 
case of a Floquet system with time-periodic Hamiltonian 
H{t) = H{t + t). This captures the non-equilibrium 
situations depicted in Fig. [T] In the case of the peri¬ 
odically driven system, we encounter the Floquet states 
|(()i(t)) = e“*®’‘/^|i(t)), which are quasistationary (i.e. 
time-periodic) solutions of the dynamics generated by 
H{t) |54H5B] . Here, |i(<)) = \i{t + r)) denotes time- 
periodic Floquet modes while Si are the quasienergies, 
which are defined modulo the energy quantum hui with 
angular driving frequency co = 2 tt jr. We start with the 
single-particle equations. In Sec. m we will then gen¬ 
eralize to the many-body case. 


A. Single-particle master equation 

We consider the time evolution of the density operator 
/3 in a single-particle system. In the weak-coupling limit, 
where the full rotating-wave approximation is valid, this 
time evolution is governed by a master equation of Lind- 
blad type |3H] i which in the interaction picture reads 

( 1 ) 

Here {A, H} = AB + BA denotes the anticommutator. 
The indices enumerate the energy eigenstates of the au¬ 
tonomous system, or the Floquet states of the period¬ 
ically driven system. In practice, we will restrict the 
number M of participating single-particle states to be 
finite. The dissipation causes transitions from eigen¬ 
state |j) to eigenstate |i) according to the jump operator 
Lij = |i)(jj, where Rij is the corresponding transition 
rate. This description is valid in the weak-coupling limit, 
where the level broadening HRij due to the transitions 
is much smaller than the (typical) energy separation of 
neighboring (quasi)energy levels in the spectrum of the 
system. The characteristic time scale rg of the unitary 
dynamics is then much smaller than the time scale tr 
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of the dissipative relaxation, rg <§; which allows to 
employ the full rotatine-wave approximation leading to 

Eq. 0 [MSI]. 

Since the resulting Lindblad equation 0 is diagonal 
in the basis of states |i), the dynamics of the occupation 
probabilities pi = {i\p\i) decouples from the off-diagonal 
elements of the density operator, which decay as one ap¬ 
proaches the steady state. The dynamics of the diagonal 
elements are described by the Pauli master equation 

Pi{t) = ^ ~ RjiPi{t)] ■ (2) 

3 

The terms of the sum correspond to the net probabil¬ 
ity flux from states j to state i. The uniqueness of the 
steady state p = '^^Pi\i){i\^ obtained by requiringpi = 0, 
is guaranteed by the Frobenius-Perron theorem, which 
holds if every state is connected with all the other states 
by a sequence of transitions with non-vanishing rates m- 

For the weak coupling to the environment considered 
here, the rates Rij in Eq. 0 can in general be determined 
in the Born-Markov (Floquet-Born-Markov) approxima¬ 
tion for autonomous (time-periodically driven) systems. 
We will consider that a bath is given by a collection of 
harmonic oscillators a with angular frequency Wq and 
annihilation operator ba, described by the bath Hamil¬ 
tonian Hb — ^ab]J}a- The bath is in thermal equi¬ 
librium with temperature T and coupled to the system 
via the Hamiltonian Hsb = v Cq(6)j -I- ba), where Ca 
are the coupling parameters and v a coupling operator 
acting in the state-space of the system. 

Within the Floquet-Born-Markov approximation, the 
rates for the driven system are given by Fermi’s golden 
rule mM, 

°° 9-77 

m——oo 

(3) 

Here Vji{m) = ^f^e 

Fourier coefficients of the coupling matrix elements, 
where the index m accounts for the absorption or emis¬ 
sion of \m\ energy quanta Huj due to the driving. The 
quantity 

= (4) 

is the bath correlation function, determined by the in¬ 
verse temperature /? = 1/T (the Boltzmann constant is 
set to one) and the spectral density 


Here now denote the matrix elements 

of the coupling operator of heat bath b with respect 
to the eigenstates |i) with energy E^. The rate is fur¬ 
ther characterized by the correlation functions gb{E) = 
Jf,(F)[exp(/3bF) — 1]“^ of both baths, with spectral den¬ 
sity Jh{E) and inverse temperature /?{,. 

Later we will see that the rate-asymmetry matrix 

Aij = Rij — Rji ( 7 ) 

plays a major role since many properties of the system de¬ 
pend on this matrix only. In the time-periodically driven 
case it reads 

OO 

ra——oo 

At'’ =AA - AA = 

( 8 ) 

whereas for the autonomous system one has 

&e{i.2} 

AfT =rTT - rT’ = X AT- Ei). (9) 

Note that the rate-asymmetry matrix is independent of 
the bath temperature(s). 

In contrast to equilibrium, a non-equilibrium steady 
state can retain a constant energy flow through the sys¬ 
tem. For the periodically driven system, the transition 
described by the rate R^jf^ causes a change of the bath 
energy by Ei — Ej -I- mhuj. The total energy flow from the 
system to the bath is thus given by 

Q{t) =^{Ei- Ej + mhu3)R!'^\i{t). (10) 

ijm 

Note that also pseudotransitions described by rates 
RTT^'^’ contribute to the heat flow [^. These transi¬ 
tions change the state of the bath, but not that of the 
system. For the autonomous system the energy flow into 
bath b reads 

Qbit)=J2iE^-E,)RlTp^it). (II) 

B. Master equation for the ideal quantum gas 


J{E) = Y, cl[6{E-nu3a)-S{E+hu3a)] = -J{-E). (5) 

CK 


We will assume Ohmic baths characterized by a spectral 
density that increases linearly with E, J{E) oc E. 

In the autonomous system, Eq. 0 simplifies to 


Rji — 


E 

bE{l,2} 


Rji 7 


2tt 

n 


RlT = '^AA'^9biE,-E,), ( 6 ) 


We now generalize the single-particle problem to a gas 
of N indistinguishable, non-interacting particles. In our 
approach we assume the total particle number N to be 
fixed, like in the canonical ensemble. For our considera¬ 
tions the canonical description poses the advantage that 
it contains the single-particle case as the natural limit 
fV = 1, and does not require to define new terms describ¬ 
ing the particle exchange with the bath. 
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The many-body Hilbert space is spanned by Fock 
states enumerated by the occupation numbers of the 
M single-particle states, n = (rii, 712 ,..., um)- To ob¬ 
tain the many-body rate equations we replace the single¬ 
particle jump operators = |i)(jj in Eq. Q by their 
Fock-space representation 

~ (12) 

Here denotes the annihilation operator of a particle, 
boson or fermion, in the single-particle mode i. Quan¬ 
tum jumps still correspond to processes transferring a 
single particle from one mode to another. The validity 
of the full rotating-wave approximation is, thus, still de¬ 
termined by the single-particle problem. Moreover, the 
total particle number N is conserved by the dynamics. 

As before, the dynamics of the many-body occupation 
probabilities pn = {n\p\n) decouple from the off-diagonal 
elements, which decay over time. The corresponding 
equations of motion are now given by (see Appendix 
for details) 

Pn{i') — ^ ) (f ~b (JTlj'jTli \_RijPnji(,i') RjiPnij^'jl ; (1^) 

b 

which is the many-body generalization of the Pauli mas¬ 
ter equation ([^. Here nji = (ni ,... ,ni — 1,...,rij -I- 
1,...) denotes the occupation numbers obtained from 
n by transferring one particle from i to j. The effec¬ 
tive transition rate depends on the quantum statistics 
via the choice of tr, with ct = 1 for bosons (reflecting 
the enhancement of transitions into occupied states) and 
CT = — 1 for fermions (reflecting the Pauli exclusion prin¬ 
ciple). The classical case of distinguishable (Boltzmann) 
particles corresponds to ct = 0 ; here the transition rates 
are independent of the occupation of the final statej^ 

For the periodically driven ideal gas the energy flow 
from the system into the bath is given by 

Q{t) + <^nj)n,pn{t) 

m n ij 

-b [{ni){t) + cr{nihj){t)]. 

Tti ij 

(14) 

Analogously, for the autonomous ideal gas the energy 
flow into bath b reads 

Qb{t) =EE('®* “ + anj)niPn{t) 

n ij 

= E(^* “ [{n^){t) + cr{ninj){t)]. (15) 

ij 


^ The bosonic master equation ( |l3[ l with ct = 1, as well as the 
corresponding mean-field equation also resemble rate equa¬ 
tions that are used to describe stochastic processes in classical 
systems, as we mention them already in the introduction. 


C. Non-equilibrium steady state 


In the following we are interested in the properties of 
the steady state of the ideal quantum ras, whose density 
operator shall simply be denoted by /ojjit is diagonal in 
the occupation number basis. 




(16) 


with pn determined by solving Eq. (13) for pn — 0. The 


uniqueness of the steady state m is inherited from the 
single-particle system, since every Eock state is connected 
to every other Eock state by a sequence of allowed single¬ 
particle transitions when this is assumed for the single¬ 
particle system. 

The steady-state expectation value of an arbitrary ob¬ 
servable 0 is denoted by 


(d) = tr(/3d). 


(17) 


Expectation values that we will consider in the following 
are the mean occupations that we denote by 


ni = {ni), (18) 

with the number operator fii = a| and the two-particle 
correlations {fiifij) or, rather, their non-trivial part 

Cij = {nihj) - fiifij = {{hi - hi){hj - hj)). (19) 

For the scenarios depicted in Fig. the steady state 
of the system will be a non-equilibrium steady state. 
This can be illustrated already on the level of the single¬ 
particle problem ([^. Let us first recapitulate the case of 
thermal equilibrium. The transitions induced by a single 
bath of inverse temperature j3 in an autonomous system 
are described by rates that obey 

Rij 

This can be inferred from Eq. for the case of a single 
bath. This condition implies that the steady state, ob¬ 
tained by solving Eq. ([^ is given by the Gibbs state with 
Pi = and Z = For this equilibrium 

state, the sum on the right-hand side of Eq. vanishes 
term by term. Thus, the net probability flux between two 
states i and j vanishes. This is the property of detailed 
balance, which is characteristic for the thermodynamic 
equilibrium. 


^ Whenever we are discussing transient behavior and time- 
dependent quantities (which happens only a few times) this will 
be indicated by writing out explicitly the time argument. For 
example, p(t) denotes the time-dependent density operator or 
(6}(t) a time-dependent expectation value. Otherwise, i.e. when 
writing p or (d), we are always referring to steady-state quanti¬ 
ties. 
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FIG. 2. (color online) Two model systems, (a) Tight-binding 
chain coupled to two heat baths of respective temperatures Ti 
and T 2 and coupling strengths 71 and 72 . (b) Tight-binding 
chain subjected to a time-periodic potential modulation at 
one end with driving strength 71 ^ and angular frequency cu 
and coupled to a heat bath of temperature T at the other end 
with coupling strength 7 . 


The rates characterizing the periodically driven sys¬ 
tem, Eq. ([^, or the autonomous system coupled to two 
heat baths of different temperature, Eq. § , are a sum 
of rates corresponding to different energy changes in the 
bath or to different bath temperatures, respectively. As 
a consequence, they do not obey the condition ( 20 ) any¬ 
more. This implies that, generally, the steady state also 
does not fulfill detailed balance anymore. While the net 
probability flux into a state z, determined by the right- 
hand-side of Eq. ([^, still has to vanish, the probability 
current from a certain state j to state i can be non-zero, 
i.e. the sum in Eq. ([^ does not vanish term by term. The 
lack of detailed balance characterizes a non-equilibrium 
steady state. In contrast to the equilibrium state, which 
is determined by the temperature of the bath only, the 
non-equilibrium steady state depends on the very details 
of the bath(s) (the temperature, the coupling operator, 
and the spectral density). This makes the computation 
of the many-body non-equilibrium steady state a difficult 
problem. However, it also offers opportunities to realize 
states with properties that are hard (or impossible) to 
achieve in equilibrium. 


D. Model systems 


ator for a particle at site i. The single-particle eigen¬ 
states |z), with z = 0,1,..., M — 1, are delocalized. They 
are described by wave functions {£\i) oc sin(fci£), with 
wave numbers ki = {i + + 1) and possess ener¬ 

gies Ei = — 2 J cos(fci) between —2J and 2 J. As sketched 
in Fig.j^a), the chain is coupled to two baths, on the left 
and right end of the chain. The left (right) bath is locally 
coupled to the first (next-to-last) site of the chain via the 
coupling operators ui = 7 ic|ci and V 2 = 72 c]v^_iCm-i, 
respectively}^ This coupling describes a bath induced 
fluctuation of the on-site energy. The steady state will 
depend on the coupling strength only through their rel¬ 
ative weight 72 / 71 , while their absolute weight deter¬ 
mines how fast the system relaxes. The temperatures of 
the baths are different from each other. We will, more¬ 
over, mainly focus on the interesting case where one of 
the baths is population inverted. For such a situation 
the notion of the single-particle ground state becomes 
meaningless, allowing for fragmented Bose condensation 
with multiple condensates [23] , see Section IV below. We 
model the population inverted bath by a negative tem¬ 
perature T 2 < 0 and a spectrum that is bounded from 
above (wq, < 0 ). 

The second model system is also given by a tight- 
binding chain of M sites. However, instead of coupling it 
to a second bath, the chain is periodically driven in time. 
Its Hamiltonian is given by 


M-l 

H{t) = -j'^ {c\ci+i + h.c.^ -b cos(a;t)c[^CM, 

i.=i 

( 22 ) 

with the dimensionless driving strength 7 ^ and angular 
frequency ui. The coupling to a bath of inverse tempera¬ 
ture P is realized via the coupling operator v = 7 c|ci, as 
depicted in Fig. [^b). The steady state will depend on 
the dimensionless driving strength 7 ,^, which determines 
the single-particle Floquet modes and the structure of 
the rate matrix Rij. However, the coupling strength to 
the heat bath 7 has no impact on the steady state, but 
rather determines how fast the system relaxes. 

Finally, as a third model, we consider a system of M 
single-particle states with the transition rates Rij given 
by uncorrelated random numbers, independently drawn 
from an exponential distribution 

P(i?,,) = A-^exp(-Ai?,,)- (23) 


Throughout this paper, we will illustrate our findings 
using three different model systems. Let us briefly define 
them here. Note that our results are not limited to these 
example systems. 

The first model system is a tight-binding chain of M 
lattice sites. It is described by the Hamiltonian 

M-l 

H = -J^ {c\ci+i -b h.c.), ( 21 ) 

wherein cn (cj) denotes the annihilation (creation) oper- 


The parameter A controls the time scale of the relaxation, 
but does not influence the steady state. The diagonal el¬ 
ements Rii can be set to 0 as they drop out of all relevant 
equations (such as Eq. (§). This choice of rates clearly 
models a non-equilibrium situation, since detailed bal¬ 
ance is violated almost surely. It is motivated by the rates 


^ We avoid the choice of coupling the second bath to the last site 
M, since for such a symmetric configuration the generic effect of 
fragmented Bose condensation m is absent. 
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computed for fully chaotic periodically driven quantum 
systems coupled to a heat bath [59] . A concrete example 
is given by the kicked rotor coupled to a bath which is 
discussed for single particles in Ref. m and for many 
particles in the supplemental material of Ref. [23] . 


III. METHODS 


In this paper we are interested in the properties of non¬ 
equilibrium steady states (161 of driven dissipative ideal 
quantum gases of N particles, described by the master 
equation © with jump operators or, equivalently, 
by the rate equation (131. Even though the particles are 
non-interacting, finding the steady state is a true many- 
body problem. Unlike in equilibrium, the many-particle 
solution cannot be obtained from the single-particle solu¬ 
tion in a straight-forward manner. This is a consequence 
of the interaction with the bath and reflected in the fact 
that the right-hand side of the master equation (§ is 
quadratic in the jump operators (121 and, thus, quar- 
tic in the bosonic or fermionic field operators As a 
consequence, equation (131 quickly becomes intractable 
when the particle number is increased. Therefore, it is 
crucial to develop and apply suitable methods for the ap¬ 
proximate treatment of the problem. This shall be done 
in this section. 

In the following, we will first describe quantum-jump- 
type Monte-Carlo simulations based on averaging over 
random walks in the classical space of sharp occupation 
numbers. This numerical method is quasi exact (the sta¬ 
tistical error is controlled) and allows for the treatment 
of moderately large systems. In order to treat even larger 
systems and to obtain an intuitive picture of the dynam¬ 
ics, we will then describe a mean-field theory, which will 
be based on a description in terms of the mean occu¬ 
pations hi. Finally, we augment the mean-field theory 
by taking into account fluctuations given by non-trivial 
two-particle correlations. 


A. Monte-Carlo simulations 

Quantum-jump Monte-Carlo simulations ISIl M are 
an efficient method for computing the time evolution of 
open quantum systems described by a Markovian mas¬ 
ter equation of Lindblad form. Instead of integrating the 
time evolution of the full density matrix, the method is 
based on integrating the time evolution of single states 
(the Monte-Carlo wave function). In doing so, the dissi¬ 
pative effect of the environment is included by interrupt¬ 
ing the continuous time evolution by a sudden quantum 
jump, described by one of the jump operators. When 
such a quantum jump occurs, and which one, is drawn 
from a suitable probability distribution. The time evolu¬ 
tion of expectation values can then be obtained by aver¬ 
aging over an ensemble of Monte-Carlo wave functions. 



0.0 0.2 0.4 0.6 0.8 ^ 1.0 


FIG. 3. (color online) Time evolution of the mean occupa¬ 
tions hi(t) for one realization of the random-rate model for 
M = 5 states and N = 100 particles. Time is measured in 
units of the inverse mean rate A [see Eq. (23l]. Initially, each 
single-particle state is occupied with the same probability. 
The thin lines are obtained from a single Monte-Carlo wave 
function, the intermediate lines from an ensemble of L = 1000 
Monte-Carlo wave functions, and the thick lines from mean- 
field theory. The mean-field results show small systematic 
deviations from the Monte-Carlo result. 


The error depends on the ensemble size and can, in prin¬ 
ciple, be made arbitrarily small. 

When treating the master equation ([^ with jump op¬ 
erators (12) we encounter a convenient situation. The 
dissipation can be described by jump operators (12) that 
transfer a particle from one single-particle eigenstate (or 
Floquet state) to another one, i.e. between two states of 
sharp occupation numbers n. At the same time, these 
occupation numbers are conserved by the evolution gen¬ 
erated by the system Hamiltonian, since we are dealing 
with a system of non-interacting particles. Therefore, the 
time evolution is exhausted by taking into account quan¬ 
tum jumps. This corresponds to a random walk in the 
classical space spanned by the Fock states \n) (not their 
superpositions). The Monte-Carlo wave function \n{t)) 
jumps between Fock states [rifc), in which it resides for 
time intervals of length tk , 


\n{t)) = [rife) with k such that Tk-i < t < Tk, (24) 


where Tk = Y!i=iii- 

We use the Gillespie algorithm [53] in order to com¬ 
pute the time evolution. At the beginning, the system is 
prepared according to the chosen initial conditions. Then 
the algorithm alternates between the following two steps, 
(i) The time interval tk determining how long the system 
will remain in the current state is drawn randomly from 
an exponential distribution P(tk) oc exp{—tk/t{nk)) with 
mean dwell time 


t(nfc) 


I 


(25) 


(ii) The new state with occupation rik+i is drawn ran¬ 
domly with branching probability reflecting the many- 
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body transition rates + fjnj)ni. Since only single¬ 

particle jumps are involved in Eq. (13), the next state is 
obtained from the current state by transferring a particle 
from a randomly drawn departure state i to the randomly 
drawn target state j. This single-particle jump has the 
probability 


state of ideal quantum gases, it is desirable to use also 
analytical methods. One of them is a mean-held descrip¬ 
tion of the system in terms of the mean occupations fii 

m- 

The time evolution of the mean occupations is given 
by the equations 


P{i jjUk) = -I- anj)ni. (26) 

These two steps are repeated until Tk = X]f=i exceeds 
the desired evolution time tgn- 

From an ensemble of L Monte-Carlo wave functions 
|n(“)(t)) labeled by a = 1,2 ,...,L, one can then com¬ 
pute the expectation value of an observable d. 


(d)ensemble(i) =\ ^ | d|n(“) (t)) . (27) 

Q —1 

Figure shows the time evolution of the mean occupa¬ 
tions {ni){t) for N = 100 particles on M = 5 states, for a 
single Monte-Carlo wave function (thin lines) and for an 
ensemble with L = 1000 (intermediate lines). One can 
clearly observe the relaxation to a steady state reached 
after a relaxation time of ~ 0.5. Slight temporal huc- 
tuations observed for times t > decrease with ensemble 
size L. The mean-held theory (thick lines) described be¬ 
low predicts the occupations rather well, but with small 
systematic deviations from the Monte-Carlo result. 

When computing steady-state expectation values (6), 
the effect of temporal fluctuations can be reduced by 
combining ensemble averaging with time averaging. 




(a) I -vi (ck) 

\o\nl 




ol— 1 k 




(28) 


Here it is useful to constrain the inner sum to A: > 
with such that T)_(c.) > tr, in order to exclude the 
transient relaxation process from the time average. Since 
we assume that every state is connected with all the other 
states by a sequence of transitions with non-vanishing 
rates, one can obtain accurate steady-state expectation 
values from a single Monte-Carlo trajectory, provided tgn 
is sufficiently large so that the system forgets its initial 
state after a certain correlation time. Averaging over a 
long time is, therefore, equivalent to averaging over an 
ensemble. We determine these uncertainties according 
to the Gelman-Rubin criterion , generally setting the 

relative uncertainties below one percent (small enough 
to make statistical fluctuations barely noticeable in any 
figure). For a bosonic system, this allows us to access 
particle numbers N ~ 10^ for M = 100 single-particle 
states. 


B. Mean-field theory 

In order to treat even larger systems and to gain some 
intuitive understanding of the non-equilibrium steady 


n^{t) =tr 


dt 


dC 


= { K (^)] 

3 ^ 

- Rji[ni{t) + a{ninj){t)\ | 


(29) 


for all i (see Appendix]^. Here we encounter the typical 
hierarchy: The time evolution of single-particle correla¬ 
tions (expectation values of operators that are quadratic 
in the field operators) are governed by two-particle corre¬ 
lations (expectation values of operators that are quartic 
in the field operators). The evolution of the latter will in 
turn be determined by three-particle correlations and so 
on. 

In order to obtain a closed set of equations in terms 
of the mean occupations, we employ the factorization 
approximation 

{hinj){t) = ni{t)nj{t) + Qj{t) « ni{t)nj{t) (30) 

for i ^ j. Here non-trivial correlations are neglected, 
Qij{t) « 0, so that two-particle correlations are approxi¬ 
mated by a product of single-particle expectation values 
as if Wick’s theorem was valid. In this way we arrive at 
the set of non-linear mean-field equations 

|i?ij%(t)[l -hcrhj(t)] 

3 ^ 

- i?jjhi(t)[l-I-(Thj (t)] |. (31) 

In the classical case of distinguishable particles, which 
can be shown to be captured by u = 0, the mean-field 
equation is exact. In this case, the equations of motion 
for the mean occupations hi(t) are of the same form as 
the single-particle master equation ([^ for the probabili¬ 
ties Therefore, in the classical system the mean oc¬ 

cupations are determined by the single-particle problem 
and read ni(t) = Pi{t)N. In contrast, for quantum gases 
of indistinguishable bosons or fermions the dynamics and 
the steady state will depend in a non-trivial way on the 
total particle number. In this case, the classical solu¬ 
tion can still be an approximate solution of the quantum 
system as long as hi ^ I for all i, so that two-particle 
correlations {fiifij) are negligible. However, as soon as 
the quantum degenerate regime is reached, where hi ^ I 
at least for some f, quantum statistics and with that the 
particle number will matter. 




The mean-field equations of motion can also be ob¬ 
tained by making a Gaussian ansatz, 


1 

Pa = ^ 


- ^ Vini 


(32) 


with partition function Z for the many-body density op¬ 
erator. For this ansatz the mean occupations are given 
by 


{ni)g = 


1 


— a 


(33) 


Thus, the M parameters defining the Gaussian state 
are determined completely by the M mean-occupations, 
rji = ln((hi)“^ + O'), as they can be obtained by solving 


the mean-held equations Eqs. (31). Non-trivial correla¬ 


tions vanish and multi-particle correlation functions can 
be decomposed into products of single-particle correla¬ 
tions determined by Wick decomposition. For the two- 
particle correlations the Gaussian ansatz gives [53] 

= OM(t+■")<*■>»+ 1) (34) 

i {ni}g{nj)g for i 7 ^ j. ^ 

for bosons (cr = 1) and fermions {a = —1). For i ^ j 
we hnd {ninj)g = {ni)g{hj)g. Therefore, starting from 


Eq. (29) and making the Gaussian ansatz for the density 


operator, we recover the mean-held equations of motion 


([31) with ni{t) = {ni)g. 

\Vith the quantities (hf)g, the Gaussian ansatz also 
determines the huctuations of the occupations fii as well 
as of the total particle number N = ^ ■ hi. One hnds 

{{hi - {hi)g)^)g = {hDg - {hi)l = {hi)g + a{hi)l (35) 
and 

i 

+ 51 o') ( 36 ) 

hi hi 

= (^*)9)^)g- (37) 

i 

The Gaussian state does not describe a system with a 
sharp particle number, so that we can only require that 
the mean particle number obeys 


{N)g = N. 


(38) 


Fluctuations of the total particle number are an im¬ 
mediate consequence of enforcing trivial correlations 
(hihj) = hihj for i ^ j (unless also the occupations 
of the individual states i are sharp so that their num¬ 
ber huctuations (hf) — hf vanish). This can be seen from 


Eq. (|36| , where we have not yet used the properties of the 

It is intuitively clear that 


Gaussian state like in Eq. 
a sharp total particle number induces non-trivial corre¬ 
lations among the occupations. If the measurement of 


the occupation hi gives a value Ui that is smaller (larger) 
than the expectation value hi, a sharp total particle num¬ 
ber implies that the number of particles in all other states 
is given by iV — and, thus, larger (smaller) than the 
original expectation value N — hi. As a consequence, the 
probability of measuring a certain value nj of the occupa¬ 
tion hj with j i will depend on the value measured 
for the occupation hi. 

The role played by fluctuations of the total particle 
number becomes less and less important in large systems. 
Namely, the variance of the total particle number (37) is 


the sum over the variances of the occupations of individ¬ 
ual modes (35), which are intensive. Thus the fluctua¬ 


tions of the total particle number grow in a subextensive 
fashion like the square root of the system size. That 
is the relative fluctuations of the total particle number 
vanish in the limit of large systems. This is the mecha¬ 
nism underlying the equivalence of the canonical and the 
grand-canonical ensemble. There is one important ex¬ 
ception, however. This is the case of Bose-Einstein con¬ 
densation, where in a bosonic system a mode i acquires 
a macroscopic occupation. If the total particle number 
is not conserved also the number fluctuations of the con¬ 
densate mode will be as large as the number of condensed 
particles; in this case the right-hand side of Eq. (35) is 
dominated by the second term. The extensive number 
fluctuations in the condensate mode will then dominate 


the sum of Eq. (37) and give rise to extensive total num¬ 


ber fluctuations, which are non-negligible in large sys¬ 
tems. This phenomenon is know as the grand-canonical 
fluctuation catastrophe [55] . 

However, one should note that the dynamics of the 
mean occupations hi (t) described by Eq. (29) do not de¬ 


pend on the occupation number fluctuations of the modes 
(the term j = i vanishes so that (hf) does not enter on 
the right-hand side). The mean-held equations of motion 
(31) can, therefore, provide a good approximation to the 


mean occupations also in systems featuring Bose con¬ 
densation (see reference [23] )■ This can be seen also in 
Fig-i where despite the fact that half of the particles 
occupy a single mode, mean-held theory accurately de¬ 
scribes both the transient and the long-time behavior of 
the mean occupations. 

The grand-canoncial ensemble of an ideal quantum gas 
in equilibrium with inverse temperature /3 and chemical 
potential /i is described by a Gaussian density opera¬ 
tor (32) with rji = (3{Ei — p). The mean occupations 
Eq. (33) follow the Bose-Einstein (Fermi-Dirac) distri¬ 
bution for (T = 1 (cr = — 1). The grand-canonical ideal 
gas is thus described exactly within the mean-held the¬ 
ory. This can be seen explicitly by plugging the Gaussian 
state pn oc Hi (solving the mean-held equa¬ 

tion) into the full many-body rate equations (13). By 
employing condition (20), which is fulhlled in an equilib¬ 
rium situation, one can see that the sum on the right- 
hand side vanishes term by term. This implies also that 
the equilibrium state obeys detailed balance as it should. 
Deviations from mean-held theory occur as a consequence 
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of two factors, (i) the assumption of a sharp total parti¬ 
cle number and (ii) the violation of the detailed-balance 
condition ([20|. 


Both factors (i) and (ii) are independent of each other, 
as can be illustrated using two examples. The canonical 
equilibrium state with sharp particle number is charac¬ 
terized by the non-Gaussian probabilities 

= I ^ exp (- l3E,n,) if Y.^n^ = N 

" 10 otherwise 

with the partition function Zjq. This state can be ob¬ 
tained by projecting the Gaussian state onto the subspace 
of sharp total particle number N. As a consequence of the 
sharp particle number, it does not solve the mean-field 
equation, as was discussed above. However, it still obeys 
detailed balance. Namely, plugging it into Eq. (131 the 


sum on the right-hand-side vanishes term by term as long 
as the condition (20) is fulfilled. On the other hand, we 


can allow the particle number to fluctuate freely, but vio¬ 
late condition (20). Then it will generally not be possible 


to find a solution of the mean-field form (32) that solves 


the many-body rate equations (13), because the number 


of independent equations exceeds the number of param¬ 
eters "qi- In the following, we are interested in the situa¬ 
tion, where a system of sharp particle number is driven 
into a steady state far away from equilibrium, so that 
both factors (i) and (ii) are present. Here, the mean-field 
theory can still provide a good approximation, as can be 
checked by comparing it to quasi-exact results obtained 
from Monte-Garlo simulations. 

Within the mean-field approximation, the heat flow 


for the autonomous system to bath 6, given by Eq. (15), 
takes the form 

= E [1 + ■ (40) 

The heat flow from the periodically driven ideal gas into 
the heat bath ( [l4| ) reads 

Q{t) = E E (^* “ ^3 4 " ni{t) [l -I- crhj(t)] 

m i,j^i 

+ EE [uiit) + a{n1)(t)\. (41) 

m i 

Here the second sum captures the heat flow related to 
pseudotransitions [see discussion below Eq. (10)]. Their 


contribution depends on (h|) and, thus, on the occu¬ 
pation number fluctuations of the modes. However, as 
discussed above, in a bosonic system of sharp total par¬ 
ticle number and where some modes feature macroscopic 
occupation, the Gaussian expectation value {nf)g = 
{fii )j2(h,),+ l] does generally not provide a good ap¬ 
proximation for the condensate mode(s). Therefore, it 
might be useful to introduce another approximation for 
(hf) in an ad hoc fashion. Another possibility is to aug¬ 
ment the mean-field theory such that it is able to treat 


systems with sharp particle number and, thus, with non¬ 
trivial two-particle correlations. Such a method will be 
presented in the following subsection. 


C. Augmented mean-fleld theory 


By construction, the mean-field theory fails to take 
into account non-trivial two-particle correlations Qj as 
they result from having a sharp total particle num¬ 
ber and from driving the system out of equilibrium, 
so that the detailed-balance condition (20) is violated. 


The effects of a fluctuating total number of particles 
can be assessed by projecting the Gaussian state onto 
the subspace of A-particle states, Pproj oc PnPqPn 
with Pm = ft =Ar This introduces non¬ 

trivial correlations, which can be obtained from {hihj) = 
tr (Pprojhihj). However, evaluating this matrix element is 
an onerous task even within efficient algorithms (see Ap- 
pendix[C|for an example), since all A^-particle Fock states 
have to be accounted for. Moreover, such an approach 
still does not include effects related to the breaking of 
detailed balance. 

In order to include the effects of non-trivial occupation 
correlations and fluctuations by analytic means, we in¬ 
troduce an augmented mean-field theory. This approach 
includes the two-point correlation functions (h^hi) into 
the hierarchy of equations of motions. In the original full 
hierarchy, the corresponding equations of motion take the 
form 

^(hfchi) =E + Aij){hkhihj) + Rkj{hihj) 

j 

+ Rij{hkhj) - {Rjk + Rji) (hkhi) 

+ Sik [Rkjinj +a{hkhj)) 

+ Rjkifik + cr(nfchj))] } 

- Rikink + a{hkhi)) - Rki{ni E a{hkhi)). 

(42) 

Here, as well as in the rest of this subsection, we suppress 
time arguments. This equation still involves the third- 
order correlations (hfeh^hj). 

The hierarchy can be closed by assuming trivial three- 
particle correlations. For that purpose we separate the 
number operators like hi = hi-\- Ci into their mean values 
hi and their fluctuations 


Q = h^- rii with {Q = 0. 
We now approximate 


(43) 


(a 6 o) = (44) 

while allowing, in contrast to mean-field theory, for non¬ 


trivial two-particle correlations Cfei = (CfeCi) [Eq. (19)]. 
Thus, the equations of motion for the mean occupations 
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are given by 


V^k^j + Ckj] 

j 

+ ^ - Rjknk ), (45) 

j 


which is equivalent to the exact equation (29). The equa¬ 
tions of motion for the non-trivial two-particle correla¬ 


tions are obtained from Eq. (421 by employing the ap¬ 
proximation (44). It is non-1 
pendix|^for details) 


inear and reads (see Ap- 


dCfcz 

dt 


j 

4 " RkjCij 4 ” RijCkj {Rjk 4 - Rji'jCki 

4 - Ck(^6ki ^ji)i.Rkj 4 - Rjk')iklk‘klj 4 - Ckj^ 

+ {5ki - Sji){Rkjnj + Rjkfik)}- (46) 


with (T = 1 for bosons are solved by a steady state char¬ 
acterized by the mean occupations 


1 


n* = 


e/3(B,-M) - 1 ’ 


(50) 


corresponding to Eq. (33) with rji = P{Ei — ^). For 
this solution the right-hand side of Eq. (31) vanishes 
term by term, indicating detailed balance. The oc¬ 
cupation numbers (50) obtained from the non-number- 


conserving mean-field theory correspond to the exact 
grand-canonical mean occupations m and provide a 
good approximation also for the canonical ensemble with 
sharp particle number N. In the latter case, the chemical 
potential has to be chosen such that 


= iV. (51) 

i 

Assuming the states of the system to be labeled such that 


The steady state values of fik and Cfcz have to be deter¬ 
mined by solving Eqs. (45) and (46) with the left-hand- 
side set to zero. 

Within the augmented mean-field theory the state is 
not only described in terms of the mean occupations hi, 
but also in terms of non-trivial two-particle correlations 
Cfci. As a consequence, we cannot only fix the mean total 
particle number to a value N by requiring 


{N) = Y,n^ = N. 


(47) 


Also the fluctuation of the total particle number can be 
fixed to a value AN 


(N^) - {Nf = J2C^J = 

This includes the choice 


(48) 


Eq < El < E 2 ^ , (52) 


meaningful positive occupation numbers correspond to 
values of the chemical potential below the ground-state 
energy, fi < Eq. The chemical potential increases either 
when P is increased at fixed N or when N is increased at 
fixed p. 

When in a system of finite extent, with discrete en¬ 
ergies Ei, the particle number N is increased at fixed 
P, the chemical potential will eventually approach the 
ground-state energy so that Eq — /i <C Ai — Eq. Once 
this happens at a characteristic particle number N* spec¬ 
ified below, the mean occupations of the excited states 
can be approximated by 


1 

- eP{Ei-Eo) _ 1 


for i > 1. 


(53) 


AiV = 0 (49) 

for a system of sharp particle number. Whereas the 
mean-field theory was found to be equivalent to a Gaus¬ 
sian ansatz for the density operator, we cannot give an 
analytical expression for the density operator correspond¬ 
ing to the augmented theory. 


Thus, for N N* the occupations of excited states be¬ 
come independent of /i (therefore also of N) and sat¬ 
urate. The occupation of the single-particle ground- 
state still depends on the chemical potential; assuming 
P{Eo — ^) <^ 1, one finds 


1 

- P{Eq - fk) 


= N 


0 ) 


(54) 


IV. IDEAL BOSE GASES AND BOSE 
SELECTION 

In this section we discuss in detail the steady state 
of non-interacting bosonic quantum gases. Let us first 
recapitulate the case of thermodynamic equilibrium. 


A. Equilibrium and Bose condensation 


Under equilibrium conditions, where the rates obey the 
condition (20), the mean-field equations of motion (31) 


with 


iVo ~ A ^ ^P(Ei-Eo) _ I ’ 

2>1 

such that p, ~ Eg — T/Nq. All particles that cannot be 
“accommodated” in the excited states will occupy the 
ground state. This is the phenomenon of Bose-Einstein 
condensation (or, strictly speaking, its finite size precur¬ 
sor). 

In a finite system Bose-Einstein condensation is a 
crossover, occurring when N becomes comparable to the 
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FIG. 4. (color online) Mean occupations versus total number of bosons for the steady state of a tight-binding chain of M = 20 
sites and tunneling parameter J > 0. The data is obtained from mean-field theory (thick solid lines), asymptotic mean- 
field theory (dashed lines), augmented mean-field theory (thin solid lines), and exact Monte-Carlo simulations (crosses), (a) 
Equilibrium situation, the chain is coupled to one bath of temperature T = IJ. (b) The chain is driven away from equilibrium 
by two heat baths of different positive temperature (Ti = IJ and T 2 = 0.5J), coupled to the first and the next to last site 
with 71 = 72 . (c) Same as in (b), but now the second bath is population inverted and described by the negative temperature 
T 2 = — J. The color code is the same as in panels (a) and (b), where the occupations decrease with increasing energy, (d) 
The chain is driven away from equilibrium by a periodic potential modulation at the last site with amplitude 7 ^^ = 2.3J and 
frequency Two = 1.5J. The Floquet states are colored like the stationary states (a-c) from which they evolve adiabatically when 
the driving is switched on (see Fig. 14 1 . 


characteristic value N*, which is directly given by the 
depletion of the condensate, 

^ ^ X! g/3{Ei-Eo) _ I' 

2>1 

In the thermodynamic limit, defined by taking parti¬ 
cle number N and volume V to infinity while holding 
the density n = N/V a.t a. constant finite value, Bose 
condensation is a sharp phase transition. At a critical 
density ric = N*/V, the occupation of the ground state 
becomes macroscopic and the ratio Nq/N, the conden¬ 
sate fraction, assumes a non-zero value. At the transition 
Eq — fi = T/Nq becomes zero. However, Bose conden¬ 
sation does not necessarily survive the thermodynamic 
limit. For a homogeneous Bose gas of spatial dimension¬ 
ality D < 2, the ratio N*/V diverges in the thermody¬ 
namic limit due to large occupations of low-energy states. 


so that no phase transition exists. In this case Bose con¬ 
densation can still be observed as a crossover in systems 
of finite size. This is illustrated in Fig. Qa), where we 
plot the mean occupations of a bosonic one-dimensional 
tight-binding chain of M = 20 sites versus the particle 
number N. In this system M plays the role of a dimen¬ 
sionless volume V so that the density is given by the 
dimensionless filling factor n = N/M. One can observe 
a sharp crossover: For N > N* the occupations of the 
excited states saturate so that newly added particles will 
all become part of the condensate in the ground state as 
described by Eqs. (53), (54) and (55). 


B. Driven-dissipative Bose gas and Bose selection 

The other panels of Fig. show the mean occupations 
hi versus N for situations where the tight-binding chain is 
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driven into a steady state far from equilibrium, either by 
coupling it to a second bath of different temperature or by 
time-periodic forcing (see section IID). In each of these 
panels, we can again identify a sharp crossover. When 
the particle number N reaches a characteristic value N *, 
many occupations saturate as in equilibrium. However, 
as a striking effect, newly added particles can now oc¬ 
cupy a whole group of states [Fig. |^c-d)], with constant 
relative occupations among these states. These selected 
states take over the role played by the condensate mode 
in equilibrium. This phenomenon has been termed Bose 
selection |23| . It turns out to be the generic behavior in 
the ultra degenerate regime of large density at fixed finite 
system size. 


As becomes apparent from Fig. we can distinguish 
two scenarios. Either a single state becomes selected. 
This includes the case of equilibrium Bose condensation 
depicted in panel (a), but also the non-equilibrium situ¬ 
ation shown in panel (b), where a Bose gas is driven out 
of equilibrium by the coupling to two heat baths of dif¬ 
ferent positive temperature. Or multiple states become 
selected as it can be seen in panel (c) and (d), corre¬ 
sponding to situations where a system is driven out of 
equilibrium by an additional population-inverted bath of 
negative temperature or by periodic forcing. As we will 
see in the following, the essential difference between both 
scenarios is that in the situations (a) and (b) the notion 
of the single-particle ground state is still meaningful. In 
panel (b) both baths favor larger occupations in states 
of lower energy and thus the largest occupation occurs in 
the ground state. This is not the case anymore for the 
situations (c) and (d). The population-inverted negative 
temperature bath of the system of panel (c) favors larger 
occupations in states of higher energy counteracting the 
effect of the positive-temperature bath. For the periodi¬ 
cally driven system of panel (d), the quasienergies of the 
single-particle Floquet states are determined modulo huj 
only, so that a ground state is not even defined. 


Within the scenario of having multiple selected states 
we can, furthermore, distinguish two possibilities. For 
that purpose we have to consider systems of a large num¬ 
ber of states M. In Fig. [^we plot the mean occupations 
for two systems with M = 100 states. Panel (a) corre¬ 
sponds to one realization of the random-rate model and 
panel (b) is obtained for a tight-binding chain coupled 
to a second population-inverted bath like in Fig. |^c). 
For the random-rate model (a) the number of selected 
states Ms is of the order of the system size M, roughly 
half of the states become selected for sufficiently large N. 
This implies that none of the selected states acquires a 
macroscopic occupation of the order of the total parti¬ 
cle number. For the tight-binding chain (b) we find that 
the number of selected states Ms is still of the order 
of one, namely three states are selected. As a conse¬ 
quence, each selected state acquires a macroscopic oc¬ 
cupation of the order of the total particle number and 
hosts a Bose condensate. This corresponds to fragmented 


10 ® 
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FIG. 5. (color online) Mean occupations versus total number 
of bosons for (a) one realization of the random-rate model 
with M = 100 states and (b) a tight-binding chain of M — 100 
sites coupled to two heat baths, namely one with temperature 
Ti = lOJ at the first site and a population inverted bath 
described by the negative temperature T 2 = — lOJ at the 
fifth-last site, with equal coupling strength, 71 = 72 [see inset 
of Fig.gc)]. 



Bose condensatioiQ which is therefore a generic situa¬ 
tions for driven Bose gas, unlike in equilibrium where 
this requires a rare ground state degeneracy. Thus, all in 
all we can distinguish three generic types of Bose selec¬ 
tion occurring in the ultradegenerate regime of driven- 
dissipative ideal Bose gases: standard Bose condensation 
where a single state acquires a macroscopic occupation, 
fragmented Bose condensation where a small number (of 
order one) of selected states each acquires macroscopic 
occupation, and the selection of a large number of states 
with non-extensive individual occupations that together 
attract most particles of the system. 

In the following we will provide a theory for Bose 
selection based on mean-field theory in the asymptotic 
limit of large N. It can be viewed as a generalization of 
the Eqs. (53), (l54| and (55) describing equilibrium Bose 


condensation to the case of driven-dissipative ideal Bose 
gases. Later on, also effects beyond mean-field will be 
discussed in terms of the augmented mean-field theory. 


C. Asymptotic mean-field theory 


A theoretical description of Bose selection can be based 
on mean-field theory, given by Eq. (31) with cr = 1. For 


^ Note that the system does not feature a single condensate 
in a state being a coherent superposition of the highly oc¬ 
cupied selected modes, but independent condensates in each 
mode. Namely, according to the Penrose-Onsager criterion Bose- 
Einstein condensation is defined by a macroscopic eigenvalue of 
the single-particle density matrix (ajuj } EH. In the situation 
discussed here, the off-diagonal elements of {ata^} are negligible 
as a consequence of the weak coupling to the bath. Therefore, 
each macroscopic mean-occupation Uj = corresponds to a 

macroscopic eigenvalue of the single-particle density matrix and 
an independent Bose condensate. 
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the steady state this equation reads 


(57) 


3 

for all i. Since Bose selection occurs in the asymptotic 
limit of large densities, it appears natural to approximate 

1 + hfe « hfc (58) 

in this equation. One then obtains the equations]^ 

0 Tij Rji^Tlj Tli ^ (^9) 


3 


3 


One can immediately see that some of the mean occu¬ 
pations fii have to vanish on this level of approximation. 
Namely, if we assume that a subset S of single-particle 
states possesses non-zero occupations, these states have 
to obey the linear equations 


0 = ^ TlijUj, i G S, 
3es 


(60) 


which directly follow from Eq. (59). However, with¬ 


out fine-tuning of the skew-symmetric asymmetry matrix 
Aij = —Aji^ these equations have a solution only if S con¬ 
tains an odd number of states (since a skew-symmetric 
matrix generically possesses an eigenvalue zero only when 
acting in an odd-dimensional space). Moreover, even if a 
formal solution can be found for a certain set S, it is not 
guaranteed that this solution will correspond to physi¬ 
cally meaningful solutions, where all occupation nnmbers 
are non-negative. Both conditions constrain the set 5, so 
that generically it will not contain all states. Those states 
contained in the (yet to be determined) set S correspond 
to the Bose selected states. 

In order to compute the occupations of the non- 
selected states, we have to include another level of ap¬ 
proximation. For that purpose we use that the occupa¬ 
tion of a non-selected state is determined predominantly 
by transitions from or into selected states. The large 
occupations of the selected states enhances the corre¬ 
sponding rates with respect to the rates for transitions 
from or into other non-selected states. Thus, neglecting 
transitions among non-selected states and still assuming 
rij + l « Uj Vj G S, from Eqs. (57) for non-selected states 
i we obtain 


5z - 1 


with gi = 


^3i''^3 


i i S. (61) 


^ It is interesting to note that these equations correspond to the 
conservative Lotka-Volterra equations fii = fii they 

are used to model population dynamics. Indeed, for fully con¬ 
nected rate matrices, the selected states correspond directly to 
those species that will not be extinct, but survive |46l 1471169| . 
Differences appear, however, for not fully connected rate matri¬ 
ces as will be discussed at the end of subsection IIV11 


This approximation is reminiscent of the Bogoliubov ap¬ 
proximation 120] for the weakly interacting Bose gas, 
where interactions among non-condensed particles are 
neglected. 

The set S has to be chosen such that physically mean¬ 
ingful occupations 


rii > 0 


(62) 


are obtained for all i, [i.e. both for the selected states, 
whose relative occupations are determined by Eq. (60), 
and for the non-selected states, with the occupations 
given by Eq. @]- We will prove in the following sub¬ 
section |1VE that there exists a unique set S for which 
condition (62) is fulfilled. Thus, the problem to be solved 
does not simply consist in solving Eqs. (60) and (61) for 
a given set S. It is rather the task of finding both the oc¬ 
cupations fii and the set S, for which the relations (60), 
(61), and (62) are fulfilled. 

By identifying the states of the set S with the selected 
states, we can now explain the major features of the re¬ 
sults presented in Fig.|^ One observation is that for large 
N the relative occupations among the selected states be¬ 
come independent of N. This is explained by the fact 
that these relative occupations are determined by the set 
of linear Eqs. (60), which does not depend on N. A 
second observation is that the occupations of the non- 
selected states saturate in the limit of large N. Such a 
behavior is predicted by Eq. (61), where the gi are deter¬ 
mined by the TV- independent relative occupations of the 
selected states. This implies also that the total occupa¬ 
tion of the selected states, 


= = (63) 

ies i^s 

grows linearly with TV. Finally, we can estimate the char¬ 
acteristic particle number TV* at which the crossover to 
Bose selection occurs to be given by the depletion of the 
selected states, i.e. by the total number of particles in 
non-selected states. 




i^S 


- 1 " 


(64) 


The set of selected states is determined completely by 
the rate-asymmetry matrix Aij. Namely this matrix de¬ 
termines not only the relative occupations among the se¬ 
lected states via Eqs. (60), but also the sign of the occu¬ 
pations of the non-selected states, which have to be 
positive. The latter can be seen by writing Eq. (61) as 




- 1 


E 




s. 


(65) 


Here the numerator is always positive, since both the 
rates Rij and the occupations fij are positive, and the 
sign of the denominator is determined by since it 
depends on the relative occupations among the selected 
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states, which are determined by Aij via Eqs. ( |60| . The 
fact that the rate-asymmetry matrix Aij, given by Eq. ([^ 
or by Eq. does not depend on the bath tempera- 
ture(s), implies that the set of selected states S also does 
not depend on the bath temperature(s). However, the 
occupations (611 of the non-selected states are tempera¬ 
ture d^endent, as Rij appears on the right-hand side of 


and (641 are generalizations of Eqs. (551 and (56), re- 


Eq. (65). This implies that both the total number of par¬ 


ticles in selected states Ns [Eq. (63 )] as well as the char¬ 
acteristic particle number N* [Eq. (64)] at which Bose 
selection sets in depend on the bath temperature(s). 

Finding the set of selected states S is generally a 
non-trivial problem. A brute-force algorithm would go 
through all possible sets containing an odd number of 
single-particle states, whose number grows exponentially 
with the number of modes M, until the desired set S is 
found. An efficient algorithm for finding S will be pre¬ 
sented in subsection |IV G| below. Already the question of 
how many states will be selected is not straightforward 
to answer, apart from the fact that (without fine tuning) 
it is always an odd number. 

A special case is the scenario of having a single se¬ 
lected state fc, corresponding to standard Bose conden¬ 
sation. Here, the occupations of the non-selected states 


spectively. However, the fact that the relative occupa¬ 
tions among the selected states and, even more, also 
the set S of selected states have to be determined adds 
an additional layer of complexity to the theory of non¬ 
equilibrium Bose selection. 


D. Systematic high-density expansion 

In this subsection we show that the asymptotic mean- 
field theory described in the previous subsection corre¬ 
sponds to the leading orders of a systematic expansion in 
the inverse total particle number N~^. This implies that 
it correctly captures the mean-field result in the limit of 
large N. 

Let us expand the mean occupations as a series in pow¬ 
ers of the inverse particle number N~^ 


,(i) 


(61) reduce to the simple expression 


Ui = Nvi + 

and require 




■ + ■ 




( 68 ) 


(69) 


Rki/Rik — 1 ’ 


i ^ k. 


( 66 ) 


The fact that these occupations must be positive reveals 
that this scenario occurs when the state k is ground-state- 
like in the sense that for all states i the rate Rki from i 
to k is always larger than the backward rate Rik, 


Rki — Af^i > 0 Vf ^ k. 


(67) 


The term “ground-state-like” refers to the situation of 
thermal equilibrium, where the relation (20) implies that 
the condition (67) is fulfilled for k being the ground 


state. These arguments reveal why we find a single se¬ 
lected state for the tight-binding chain which is driven 
between two heat baths of different positive temperature 
[Fig. I^b)]. In this situation the notion of the single¬ 
particle ground state still remains meaningful even away 
from equilibrium. This is generally different when the 
system is coupled to a population-inverted bath described 
by a negative temperature, like in Fig. ®c), or in a pe¬ 
riodically driven system, like in Fig. |^d). In the former 
case the condition (67) cannot be expected to hold for k 


being the ground state and in the latter case the ground 
state is not even defined (since quasienergies are deter¬ 
mined modulo huj only). 

We can compare our theory to the theory of equilib¬ 
rium Bose condensation as it was reviewed in subsection 
|IV A[ First of all, we would like to note that the equilib¬ 
rium situation is contained in our asymptotic mean-field 
theory as a special case. Namely, the equilibrium expres¬ 
sion (53) for the excited-state occupations is reproduced, 
when the relation (20) is plugged into Eq. ([6^. Gener¬ 


ally, our Eq. (61) generalizes Eq. (53); likewise Eqs. (63) 


for the leading order as well as for the corrections of or¬ 
der r > 1. These requirements ensure that the mean 
total particle number is given by when the series is 
truncated after a certain order r. Such an expansion is 
equivalent to an expansion in the inverse pa rticl e density 


n ^ = M/N. We can now plug the ansatz (68) into the 
mean-field Eqs. ((5~ 


0 A,jVj 

3 

1 

N ^ 

3 

1 


-)- ^ Rjik'i 




( 1 ) 




-b 


7V2 


E[ 


Nijl3j Rjil3, 


( 1 ) 


3 

-\-A. 




( 2 ) I 

Vi Vi + V, 




+ 


( 2 ) 


+ 0 


1 

iW 


(70) 


and ask that all terms that correspond to the same power 
of N vanish independently. In this way we get a hierarchy 
of equations determining the coefficients of the expansion 
(68) order by order. 


Collecting the terms of the leading order gives rise to 
a set of equations for the leading coefficients Vi. These 
equations take the form of Eqs. (59), but with hi replaced 
by V,, 


q = u,y a 


(71) 
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Repeating the arguments of the previous section we see 
that the leading-order coefficient is non-zero only for a 
(yet to be determined) set of selected states S, so that 

Iy^=0, iiS, (72) 


and 


0 — Aij Vj , 

j&S 


i G S. 


(73) 


The next order determines the coefficients i/) 


( 1 ) 


Thanks to Eq. (72) the coefficients of the non-selected 


states are not coupled to each other and depend on the 
leading-order occupations of the selected states only, so 
that we arrive at the simple expression 


,(i) 


V- = — 


^jeS 


5. 


(74) 


This expression directly corresponds to Eq. (61), but with 
hi replaced by t'i for the selected and by the non- 

selected states. The leading corrections to the occupa¬ 
tions of the selected states appear in the same order and 
can be determined by solving the linear equations 


j&S 


,(i) 






“t” ^ ^ 


3tS 


( 1 ) 


13''j 


), i^S, (75) 


where we used 0 = [^dS- 0 and ([T^], and 

taking Eq. (691 for r = 1 into account. Higher orders in 
the expansion (70) can become relevant when some rates 
vanish as discussed in Sec. lYHl 

Truncating the 1/iV expansion after the first order, one 
obtains 


.,( 1 ) 


( 1 ) 


for i G 5 
for i ^ S. 


(76) 


However, asymptotically in the limit of large N, it will be 
sufficient to take into account only the leading contribu¬ 
tions, so that the mean occupations can be approximated 
as 


ViN for z G 5 
for i ^ S. 


(77) 


This corresponds to the approximation of the previous 
subsection, apart from the slight difference that, previ¬ 
ously, we normalized the total occupation of the selected 
states Ns to the first-order result Ns^'^ = + 

= N — Yl,i^s implicit in Eq. (63) and 

corresponds to the approximation 


for z G 5 
for z ^ 5, 


(78) 


This normalization, which for finite N takes care of the 
fact that the leading contributions to the occupations of 
the selected and the non-selected states stem from differ¬ 


ent orders, is thus a compromise between Eq. (76) and 
Eq. ( |77| ). For large but finite N it is better than Eq. (77), 
since it produces the correct total particle number, but 
it does not require to compute corrections for the 


selected states that enter Eq. (76). Therefore, we will 


use Eq. (78), corresponding to the asymptotic theory as 


it was presented in the previous subsection, in the follow¬ 
ing. In the asymptotic limit fV —>■ oo all three expressions 
([T^, and ([T^ are, of course, equivalent. 


The requirement of having a positive particle number 
in the asymptotic limit of large N is given by 


r-j > 0 

,(i) 


for z G 5 
' > 0 for i ^ S. 


(79) 


In order to find a compact formulation of finding an 
asymptotic solution obeying this condition it is conve¬ 
nient to introduce the numbers ^ ^b^i- Accord¬ 


ing to Eq. (73) they vanish for z G 5, while Eq. (74) 


tells us that they should be negative to ensure positive 
occupations of the non-selected states. The problem of 
finding an asymptotic mean-field solution can, therefore, 
be reduced to the problem of finding a set S of selected 
states and numbers Ui and /Zi such that [23] 


/i* = 


y Aij Vj with 


Vi> Q and /z^ = 0 for z G S, 
Vi = Q and /z^ < 0 for z ^ S. 


(80) 

The non-generic situation with Vi = fii = 0 for some z 
corresponds to transitions, which we discuss in the next 
subsection. Before we prove that a unique set S obeying 
the relations (80) exists, let us point out that these re¬ 


lations are valid only in the case of fully connected rate 
matrices. If we allow for zero rates Rij = 0, the set of 


selected states is not determined by the conditions (80) 
anymore, as we discuss in subsection |IV I| below. 


It is interesting to note that the conditions (80) that 


determine the selected states are equivalent to those de¬ 
termining the surviving species under the dynamics of 
the Lotka-Volterra equations given in footnote El 
Differences appear for non fully connected rate matrices 
(see discussion at the end of Sec. IVI). 


E. Existence and uniqueness of the set of selected 
states 

In this subsection we provide a proof for the unique¬ 
ness and the existence of the set of selected states for 
fully connected rate matrices (which we repeat for com¬ 
pleteness from the supplemental material of Ref. [23].) In 
the following we will use the vector and matrix notation, 
with u and denoting the vectors with elements iZj and 
/Zi, respectively, and R and A denoting the rate matrix 
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and the rate-asymmetry matrix with elements Rij and 
Aij^ respectively. Let us, furthermore, decompose A like 




ASS 

A55 

A^ 


(81) 


wherein the submatrix = {Aij}ij^s denotes the rate- 
asymmetries among selected states, = 

{^ij}iiS,jGS rate-asymmetries among non-selected 
and selected states, and the rate- 

asymmetries among non-selected states. The conditions 
(80) with i G S require us to determine S such that 
A* has a vanishing eigenvalue. Skew-symmetric matri¬ 
ces generically have a vanishing eigenvalue only if their 
dimension is odd. As the square submatrix A^ of A is 
still skew-symmetric, we can immediately conclude that 
the number Ms of Bose selected states is odd. The con¬ 
ditions (80) stipulate, furthermore, that the correspond¬ 


ing eigenvector Vi, i G S has positive components. Fi¬ 
nally, the conditions for i ^ 5 tell us that this eigenvector 
should result in a vector with non-positive components 
when it is multiplied with the submatrix A^^. 

We now prove the uniqueness of the set S. Assume 
first that there exist two different sets iSi and ^ 2 , both 
leading to physical solutions and z ^2 with /xi = Ai^i 
and /X 2 = Av >2 obeying Eq. (80). Using 


1/2 fii = 1^2 ^^1 = (^2 ^^1)^ = ^1 ^2 = Ai/-^ 


= M2, 


it then follows from Eq. (80) that 

0 > V>2 f^2 > 0 . 


(82) 


(83) 


This requires that both i/J fii = 0 and v]" /X 2 = 0, such 
that S 2 C 5i and Si C S 2 , leading us to conclude that 
Si = S 2 = S. Given the set S, the homogeneous linear 
system for v generically has a single solution only. There¬ 
fore, the solution to the generic steady-state problem has 
to be unique. 

In order to prove the existence of the set S, we 
now restrict S to sets comprising an odd number Ms 
of states, according to the generic conditions described 
above. Each choice of S gives rise to a (possibly non¬ 
physical) solution 1/5 with fis = A 1 / 5 . The vector of 
signs (T with 


Ui = sign(i^i) if i e 5, 
ai = -sign(/ij) Hi ^ S, 


(84) 


distinguishes physical solutions (di = 1 for all i) from 
non-physical solutions. Here, we fix an overall sign due to 
the orientation of the vector 1^5 by the convention cti = 1 . 
Now we observe: (i) Cycling through all odd-numbered 
subsets S, each possible vector cr occurs at most once. 
Namely, if Si and S 2 gave rise to the same vector cr then 
the modified rate imbalance matrix = aiAijaj had 
two physical solutions with different selected sets iSi and 



FIG. 6. (color online) Mean occupations in response to the 
variation of a dimensionless parameter p, for a small system of 
N = 10® bosons on M = 7 states for the random-rate model. 
The rate matrix R{p) is a superpostion of two independently 
drawn rate matrices and R^^\ with the relative weight 
controlled by p, R{p) = (1 — p)R^^^ -\- pR^^^■ The results are 
obtained using mean-field theory (dotted lines), asymptotic 
theory (solid lines for selected states and dashed lines for non- 
selected states). Each color refers to a specific state. At 
each transition two states are exchanged between the sets of 
selected and non selected states. 


^ 2 , in contradiction to the previously established unique¬ 
ness of the solutions, (ii) The number of possible 

vectors cr equals the number X)ms=i 3 (ms) “ 2^^“^ 
of possible sets S. Therefore, each vector cr occurs once. 
In particular, this includes the vector with Ui = 1 for all 
i, leading to the solution with positive macroscopic and 
microscopic occupations. This guarantees the existence 
of a physical solution. 


F. Transitions 

In this section we will discuss transitions, where the 
set of selected states S changes in response to the vari¬ 
ation of a parameter p. Examples for such transitions 
can be observed in Eig. This figure shows the mean 
occupations versus the parameter p for a model de¬ 
fined by the superposition of two random rate matrices 
and R^'^\ with the relative weight controlled by p, 
R{p) = (1 — HpR^-^K One can see that in a transi¬ 

tion two states are exchanged between the set of selected 
states and the set of non-selected states, such that the 
number of selected states is odd before and after the tran¬ 
sition. Approaching a transition from the left, the tran¬ 
sition is found to be triggered by a state i"'. This state 
can either be a selected state whose occupation drops 
until it becomes non-selected at the transition (case I) 
or a non-selected state whose occupation increases until 
it becomes selected at the transition (case II). Further¬ 
more, one can observe that at the transition a second 
state becomes involved abruptly that changes from 
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log 

A 


log 

B 


FIG. 7. (color online) Four generic types of transitions, where 
the set of selected states changes from 5 = to 5 = 5^ 
when a parameter p reaches a critical value p*. In each tran¬ 
sition two states are exchanged between the sets of selected 
and non-selected states, so that the number Ms of selected 
states remains odd. When approaching the transition from 
the left, it is triggered by a state either a selected state 
whose occupation drops until it becomes non-selected at the 
transition (case I) or a non-selected state whose occupation 
increases until it becomes selected at the transition (case II). 
A second state becomes involved abruptly at the transi¬ 
tion that changes from selected to non-selected (case A) or 
vice versa (case B). This state plays the role of the triggering 
state when the transition is approached from the right. Types 
(I,A) and (II,B) form one class, since they transform into each 
other when the transition is passed in opposite direction. 





the selected to non-selected (case A) or vice versa (case 
B). When approaching the transition from the right, the 
states i"' and change their role, so that the former 
partner state v" plays the role of the triggering state. 

The four combinations of cases I or II and A or B de¬ 
fine four generic types of transitions that are depicted in 
Fig. Type (I,A) and type (II,B), where the number 
Ms of selected states is lowered or raised by two, respec¬ 
tively, transform into each other when the transition is 
passed in opposite direction. Therefore, they form one 
class. In type (II,A) transitions, which are triggered by 
non-selected states from both side, and type (I,B) transi¬ 
tions, which are triggered from selected states from both 
sides, the number Ms of selected states does not change. 
They define two distinct classes, since they cannot be 
transformed into each other. 

These observations based on Fig. turn out to be 
generic. In the following we will describe them within 
the asymptotic mean-field theory. We have already de¬ 
fined the left triggering state and its partner state, the 
right triggering state . Let, moreover, p* be the criti¬ 
cal parameter at which the transition occurs and and 
be the sets of selected state on the left-hand and the 
right-hand side of the transition, respectively [Fig. [^. 
Within the asymptotic theory a transition must occur 
when the occupation hi of a state i would change its sign 


at a critical parameter p = p*. This state i plays the 
role of the triggering state If is a selected state 
(before the transition), the transition occurs when ^'i< 
drops to zero, so that in zeroth order the occupation of 
this state becomes zero. In that case the state can, 
thus, be viewed as a non-selected state at the transition. 
If is a non-selected state, the transition occurs when 
/ii< becomes zero, so that in first order the occupation 
of this state diverges. In that case the state must, 
therefore, be viewed as a selected state at the transition. 
Thus, at the transition p = p*, which corresponds to a 
fine-tuned situation, the set of selected states contains an 
even number of states and is given by 


5 * 


5<U{f<} if i<^S< 
5<\{I<} ifi<G5<. 


(85) 


As the number of Bose-selected states has to become odd 
after the transition, one further state has to be in¬ 
volved. The set S* can also be expressed in terms of this 
partner state. 


5 * = 


5 >\{*>} 


if i> i S> 
if G 5>. 


( 86 ) 


In the following we will describe how to determine this 
partner state in order to find the set of selected states 
on the other side of the transition. 

The intricate details of the transition are encoded in 
the truncated matrix , obtained from A* = A{p*) 
by removing all rows and columns corresponding to non- 
selected states i ^ S* like in Eq. (81). According to the 
transition criteria, this matrix has at least one vanish¬ 
ing eigenvalue. As the matrix is even-dimensional and 
skew-symmetric, its eigenvalues are imaginary and come 
in pairs of opposite sign. Thus, one eigenvalue of zero im¬ 
plies another one, so that generically the kernel of A'^ 
will be two-dimensional at the transition. One vector 
lying in the kernel of A^ is given by the limiting occu¬ 
pations z/i as one approaches p* from below. We denote 
this vector by (note that this is now truncated to 
the states of S*). Analogously, there is a second vector 
from the limiting occupations as one approaches p* 
from the right, which also lies in the kernel. We will now 
establish a relation between both vectors and zz^. 

For that purpose, we introduce an interpolating vector 
zz(a) = aiz< -I- (1 — a)iz', where iz' is the element of the 
kernel of A* which is orthogonal to iz< while a is an inter¬ 
polation parameter. The occupations of the non-selected 
states (and their sign) is determined by the vector /x given 
by Eq. ( [80| . For the two possible solutions zz< and zz>, 
this vector reads /x< = ^ iz< and /x^ = A*^ ^ zz>, 

respectively. Both vectors are connected by the inter¬ 
polation /x(a) = a/x< -I- (1 — a)/x' with /x' = A'^ ^ iz'. 
Herein A*^ is obtained from A* = A{p*) as described 


by Eq. (81). Due to the selection criterion Eq. (80), 
we require physical solutions izi(a) > 0 Vx G 5* and 
Mi (a) < 0 Vx ^ 5*. Choosing the orientation of iz' conve¬ 
niently, this is fulfilled for the finite interval 0 < a < . 
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The extremal point is determined by ramping up a 
until either an element of i>'(a) or /x(a) becomes zero. 
The index of this element corresponds to the state 
and the extremal point determines the solution 
via ■ 


Exactly at the transition, the interval 0 < a < 
corresponds to physically meaningful solutions with pos¬ 
itive occupation numbers. Its extremal points describe 
the solutions and found when approaching the 
transition from the left and right hand side, respectively. 
The narrower the interval, i.e. the smaller Aa = a^, the 
more similar will both solutions and be. That 
means the smaller will be the discontinuous changes in 
the occupations of the states i ^ that are not 

directly involved in the transition, as they are visible 
also in Fig. The width Aa associated with a typical 
transition must, moreover, be expected to shrink with 
the system size. Namely, each of the M single-particle 
states of the system provides a constraint that poten¬ 
tially limits this interval, since the number of conditions 
Eq. (80) proliferates with M. So in large systems one 
cannot only expect more transitions to occur when a pa¬ 
rameter is varied, but also that the discontinuous jumps, 
which the non-participating occupations undergo at each 
transition, become smaller. 


Before moving on, let us briefly discuss the case of 
finite particle numbers N, where the sharp transition be¬ 
comes a crossover of finite width. This can be observed 
in Fig. Here we plot the mean occupations versus the 
total particle number N for a system described by the 
same rate matrix R{p) used in Fig. The five panels 
of Fig. 1^ are obtained for parameters p close to (or at) 
the transition labeled (I,B) in Fig. with the critical 
parameter denoted by p*. The first panel corresponds 
to a parameter well on the left-hand side of the transi¬ 
tion. Here asymptotically three states become selected. 
When coming closer to the transition, but still staying 
on its left-hand side (second panel), we can observe that 
a preasymptotic regime appears. Namely, at large, but 
finite N the system approaches a state with two selected 
states, before eventually in the asymptotic limit A —>■ oo, 
a third state becomes selected as well. This third state 
corresponds to the triggering state The two states 
that appear to be selected in this preasymptotic regime 
correspond to those two states that are selected at the 
transition (middle panel). The fourth panel corresponds 
to a parameter, where the transition has just been passed. 
Here (roughly) the same preasymptotic state is found, 
before asymptotically for A —>■ oo a third state joins the 
group of selected states. Now the third state is given 
by . The fifth panel is, finally, obtained for a param¬ 
eter well on the right-hand side of the transition. Here 
again no preasymptotic regime is found. The emergence 
of a preasymptotic regime close to the transition implies 
that the fine-tuned rate matrix R{p*), which gives rise 
to two selected states, provides an accurate description 
of the system within a finite interval of parameters near 
the transition. 


G. Efficient algorithm for finding the selected 
states 


In principle, finding the unique set S of Bose-selected 
states requires to sample all possible subsets, whose num¬ 
ber grows exponentially with M, until one succeeds to 
satisfy the conditions (80). Testing all sets by brute force 


quickly becomes unpractical already for moderately large 
values of M. While the mean-field occupations and espe¬ 
cially their dependence on the total particle number can 
provide some guidance, this method also quickly reaches 
its limits when M is further increased. Here we describe 
an efficient algorithm for finding the set of selected states. 
It uses the theory of transitions that we presented in the 
previous subsection. 

In order to solve the problem of finding the set of se¬ 
lected states for a given rate-imbalance matrix A, we con¬ 
struct the auxiliary rate-imbalance matrix 


^ij (p) — klij ~\~ pBi^ 


(87) 


by adding the real-valued skew-symmetric matrix H, 
weighted with the real parameter p, to the original one. 
The problem defined by the new matrix A{p) will be 
solved by a set S{p) of selected states. The matrix B is 
constructed as follows: It shall possess a cross-like struc¬ 
ture, with non-zero elements only in the column and the 
line labeled by k, 


Bij ^ik^j ^kj^i: (^^) 


SO that 


~ ^kj + bj > 0, Vj ^ k. 


(89) 


This condition, which corresponds to the relation 
ensures that for p = 1 only the state k will be selected, 
iS(l) = {k}. Relation (89) can be achieved with minimal 
effort by setting 


b^ = 


+ I + £i > 0 if Aj~i < 0 
0 otherwise 


(90) 


with arbitrary Si > 0. Our strategy will now consist 
in ramping the parameter p down from p = 1, where 
the solution 5(1) = {k} is known by construction, to 
p — 0, where we would like to know the solution 5(0) = 
5. During this ramp, we will monitor all transitions, 
i.e. changes of the set 5(p), that are happening, so that 
at the end we will arrive at the desired solution. For 
that purpose it seems favorable (though not necessarily 
required) to choose the state k such that a minimum of 
the elements bj defined like (90) has to be non-zero, and 
to choose the Si different from each other, Si ^ Sj for 
i j, m order to separate the transitions when varying 


P- 

In order to follow the state of the system during the 
parameter ramp, we take advantage of the specific way 
the matrix A{p) depends on the parameter p. Namely, 
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FIG. 8. (color online) Occupations numbers versus total particle number N close to a transition. The system is described by 
the same rate matrix R{p) as that of Fig. The parameters p used in the different panels are chosen to be close (or at) the 
transition labeled (I,B) in Fig. with p* = 0.303179 denoting the corresponding critical parameter. 


the cross structure (88) of the matrix B implies that the 
occupations of the system change in a linear fashion un¬ 
less a transition occurs: If the vector i>{po) solves the 
problem (80) for A{pq), then one has 


C>{p) = C{p) i>{po) + i>'{po){p - po) for Pa < P < Pb, 

(91) 

with a global normalization factor C{p) > 0 such that 
J2j^S{p) ^(p) ~ 1- Here the limits Pa and pb are given by 
those values of p, where the set of selected state changes 
away from S{po) in a transition. The proof of this state¬ 
ment is rather technical and delegated to appendix [E} 
where we also describe how to obtain v>'{pq). Expression 
(91) can be employed to predict the positions Pa and pb 
of the transitions as those points, where either an ele¬ 
ment Vi(p) of v(jp) or an element pii{p) of the associated 
vector fl{p) = A{p)u{p) would change sign. The label 
i of this state corresponds to the state that triggers the 
transition. 

With these ingredients, our algorithm works as follows: 
Start from p = 1, where 5(1) = {fc}, and evaluate where 
the next transition occurs when p is lowered and by which 
state it will be triggered. Next, employ the theory of 
transitions described in the previous subsection to deter¬ 
mine the partner state which at the transition also 
changes between the sets of selected and non-selected 
states. In this way the new set of selected states solv¬ 
ing A(p) after the transition has been found. Then com¬ 
pute where the next transition occurs when p is lowered 
further, iterating this procedure until p = 0 is reached. 
The time needed to find the set of selected states in this 
way scales polynomial with the system size M. For the 
random-rate model, which constitutes a rather difficult 
problem since on average half of the states are selected 
[23], we find this time to scale as ^ M“ with a « 4. This 
allows us to find the set of selected states for systems of 
up to M = 1000 states. An alternative algorithm for 
solving the problem (80) has recently been presented in 



FIG. 9. (color online) Effect of small rates for a minimal 
three-state model with rate matrix ( |94[ ). Occupations hi ver¬ 
sus total particle number N obtained using mean-field theory 
(solid lines) and asympt otic theory (dashed lines) for the rate 
matrix R given by Eq. (941, which is visualized in the inset 
(line widths reflect rates). Fnrthermore, the dotted lines show 
occupations hi obtained by the asymptotic theory for the ap¬ 
proximate rate matrix given by Eq. (95 I, where the small 
rates have been neglected. Blue, green, and red lines describe 
hi, h2, and ha, respectively. Near A ~ 10 the system ap¬ 
proaches a preasymptotic state with a single selected state, 
described by i?“, before above N ~ 10° the trne asymptotic 
state is reached, where all states are selected. 


Ref. [47] and is based on linear programming. 


H. Small rates and preasymptotic regime 

So far we have assumed strictly positive rates, Rij > 0, 
within the asymptotic theory. This assumption is reason¬ 
able in the sense that exactly vanishing rates, Rij = 0, 
can be viewed as a fine-tuned situation. However, obvi¬ 
ously, we can encounter situations where some rates are 
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much smaller than others, e.g. 


Rij — 


0{r) for (i,j)eG' 
0{er) else, 


(92) 


with G denoting the subset of pairs {i,j) with large rates 
of order r and e <C 1 quantifying the suppression of small 
rates of order er. Such rate matrices can result from a sit¬ 
uation where some modes are coupled much more weakly 
to the environment than others. Having such a situation 
in mind, in the following discussion we will consider a rate 
Rij to be small only when also its backward rate Rji is 
small too, so that also the corresponding rate-asymmetry 
\Aij\ is small 

Having some rates much smaller than others, it ap¬ 
pears reasonable to neglect the small rates in an approx¬ 
imation, 


R, 


Tja _ 

— 


Rij ioi{i,j)eG 
0 else. 


(93) 


As we will argue below, such an approximation will de¬ 
scribe the system accurately, provided the total parti¬ 
cle number N remains below a threshold A^thr associated 
with the approximation. Thus, when increasing the par¬ 
ticle number N, one might encounter the following sce¬ 
nario: First a preasymptotic state is approached, where 
the occupations are well described by the asymptotic the¬ 
ory based on the approximate rate matrix i?“-, before 
eventually the true asymptotic state of the full rate ma¬ 
trix R is reached above the threshold. This scenario can 
be observed in Fig. where we plot the occupations of 
a minimal three-state model versus N. In this model the 
rates are given by 


R = r 


=r 


{^ 

1 

2e\ 



2 

0 


e = I0-3, 

(94) 

(le 

4 

oy 



{ ° 

I 

o\ 



2 

0 

2 ■ 


(95) 


0 4 0 


This behavior resembles the preasymptotic behavior 
found near transitions that we discussed at the end of 
subsection |IVF In both cases the preasymptotic state is 
described by a fine-tuned rate matrix, either character¬ 
ized by the critical parameter or by setting several matrix 
element to zero. However, since setting several matrix el¬ 
ements to zero corresponds to the fine tuning of several 
parameters, the set of selected states of i?“ can be quite 
different from that of R. 


The appearance of a preasymptotic regime described 
by the approximate rate matrix (951 at intermediate par¬ 
ticle numbers N, as it is visible in Fig. roughly for 
10 < AT < lO^, can be explained as follows. When apply¬ 
ing the asymptotic theory. Sec. IV C[ to the approximate 
rate matrix, where small rates are neglected, we find that 
the selected state 3 acquires an occupation ~ V, while 
the occupations of the non-selected states are ^ 1. Gen¬ 
erally, the fact that the selected state(s) possesses an oc¬ 
cupation much larger than the non-selected states justi¬ 
fies the 1/N expansion (68), which underlies the asymp¬ 
totic theory. This explains why the preasymptotic regime 
is reached near N ~ 10, when V ^ 1. However, as soon 
as the factor N between the occupations of the selected 
and the non-selected states becomes comparable to the 
inverse suppression factor ^ 10^, the weak rates start 
to spoil the hierarchy of the 1 /N expansion based on the 
selected state of i?“. Namely the product of a small rate 
with the occupation of a selected state ~ reN, which was 
neglected so far, can become comparable to the product 
of a large rate with the occupation of a non-selected state 
~ r, which has been taken into account. This explains 
why for N > Vthr = 10^ the system starts to devi¬ 

ate from the solution of the approximate rate matrix (de¬ 
scribing the preasymptotic state), to approach the true 
asymptotic state determined by the full rate matrix. 

Note that allowing for zero rates, i.e. rate matrices that 
are not fully connected like i?“, can have several conse¬ 
quences for the asymptotic theory. These are discussed 
in the following subsection. 


I. Zero rates: Not fully connected rate matrices 


So far we have assumed fully connected rate matrices 
within our asymptotic theory. What happens if we allow 
some rates to become zero? This question emerges, e.g., 
when computing the asymptotic state of an approximate 
rate matrix [Eq. (93)]. First of all, in case the rate 
matrix is disconnected, so that it is not possible any¬ 
more to reach every state i from every other state j in 
a sequence of quantum jumps (and vice versa), then the 
steady state of the system is not unique anymore [57] and 
will depend on the initial conditions^ We will exclude 
this scenario from the following discussion and focus on 
situations where the rate matrix is solved by a unique 
steady state. 

In order to discuss the impact of zero rates, let us 
briefly recapitulate the situation where all states are cou¬ 
pled to all other states. In this case the coefficients Vi and 

(r) 

of the XjN expansion (68) are obtained as follows. 
First the leading coefficients Vi^ and with that the set S 


® There can also be small rates without small backward rates, e.g. 
between states with a large energy separation. Not considering 
those rates as small in the below analysis (i.e. not exploiting the 
fact that they are small) does not spoil its validity. 


^ If, by taking into account neglected rates of order er, the matrix 
is connected again, then for times longer than l/(£r) the non¬ 
unique steady states associated with will eventually relax to 
the unique steady state of the full rate matrix. 
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of selected states, have to be determined by solving the 

1 “^ fr) 

problem (80). Then the sub-leading coefficients vl can 
be obtained iteratively from the hierarchy of equations 
that results from Eq. (70) by requiring the terms of each 


power of N to vanish separately. If we denote the terms 
oc N~'^ on the right-hand side of Eq. (70) by , then 
this hierarchy of equations reads 






(96) 


for all i and for r = 0,1,2,... and with denoting 
the vector of coefficients Now the are obtained 
by solving the set of linear equations with 

the already determined Ui treated as parameters. Then 
the ^ are obtained from the set of linear equations 
= 0, with the already determined coef¬ 
ficients Vi and v\^'^ entering as parameters, and so on. 

This procedure has to be modified for non fully con¬ 
nected rate matrices. In the following discussion we will 
assume that Rij = 0 implies Rji = 0 and, thus, also 
Aij = 0, this is analogous to our assumption about 
the occurrence of small rates in the previous section. 
Let us start with the zeroth-order equation, = 

— 0- As before, we conclude that the leading 
coefficients Vi are non-zero only for a group of selected 
states i £ S, 


by Eq. (97). Thus, one has K — I parameters vs^ that 


yet have to be determined from the equations of higher 
order. 

In order to investigate the first-order equations = 
0, [see Eq. (70)] it is useful to define two groups of non- 
selected states, 


5 = 5' U S" 


( 102 ) 


such that states that are directly coupled to selected 
states via non-zero rates form the set S' and states that 
are not coupled directly to any selected state form the 
set S". For i G S' they lead to the familiar result 


vP=- 


^jes 


E 


jes 


i G 5'. 


(103) 


Note that for states i G S' that are coupled to selected 
states belonging to two sub sets Sg and Sjs (or more), 
the right-hand side of Eq. (103) depends on the ratio 


vs^lvsp, which is not determined yet. In that case the 
ratio vg^fvg^ can be obtained from Eqs. (105) below. 

The coefficients v^^'^ with i G S" drop out of the first- 
order equations (/|^^ = 0 is fulfilled trivially) and must 
be determined from the second-order Eqs. ( W below. 
The first-order equations for th e sele cted states z G 5a of 
a subset Sa simplify [with Eq. (101)] to 


je-s 

i G (S, 

(97) 

0 — ^ ^ E? Rji^i 


Vi = 0, 

* ^ 5, 

(98) 


5 

z G 5a. (104) 


The set S of selected states has still to be determined 
from the requirement that the asymptotic occupations of 
both the selected and the non-selected states are positive. 
It can consist of K uncoupled subsets Sa, 


5 = 5i U 52 U • • • U 5 k, 


(99) 


with 


Rij = 0 


for 


i G 5a, J G Sp,a ^ /?. (100) 


such that each subset Sa fulfills Eqs. (97) individually, 

= 0 , WiGSa- ( 101 ) 


E 


Aij Vj 


Without fine-tuning, a solution of '}2j(^s Aij^j = 0 is 
guaranteed as long as the number of states in each of the 
subsets Sa is odd. However, the total number of selected 
states Ms can now also be even. It is even (odd), if the 
number K of uncoupled subsets Sa is even (odd). In the 
case of fully connected rate matrices, the coefficients Vi 
were determined uniquely by the set 5, Eqs. (|9^ and 
(98), as well as by the normalization condition (|69|). For 


K > 1 this is not the case anymore. Here the relative 
occupation of a subset 5a, defined by vs^ = Ezgs 
is not fixed, since vs^jvs^ with a ^ /3 is not determined 


jeS' 


These equations determine the coefficients v^^'^ of the se¬ 
lected states i. 

Further information can be obtained by summing 
Eqs. ([^ over all states i G Sa- This gives 0 = 


EjgS' E*g5c Here, all non-selected 

states j G S' that couple only to selected states of the 
subset 5a do not contribute to the sum, since accord- 

,(i) 


ing to Eq. ( |103| ) their occupations are given by Vj ' = 
“(Ejg5c ^hus, we obtain 


0 — E! {jRji A AjiZzj , 

j^Sa-G 


Va 


(105) 


where Sa+ denotes the set of non-selected states that 
couple to the subset Sa and at least to one more selected 
state of a different sub set S g with fi ^ a. If this set 
5a+ is not empty, Eq. (105) can be used to determine 
missing relative occupations vs^ / vg^. We will argue be¬ 
low that in fact all subsets of selected states must form 
a connected cluster, where two subsets Sa and Sp are 
defined to be connected if they are coupled directly (via 
a single quantum jump of non-zero rate) to the same 
non-selected state(s). This guarantees that all relative 
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occu patio ns can be determined from Eqs. (105) 

and (103), so that the Ui can be determined completely. 


( 2 ) 

From the second-order equations ^ = 0, we obtain 


o = E 




■A,, 


,( 2 ) 




l/j + V- 'v] ' + ViV, 


(2)^ 


Vi. 

(106) 


These equations determine all the coefficients that 
have not been obtained yet, since all are coupled to 
each other (at least indirectly). For the missing coeffi¬ 
cients of states i G S" they simplify further to 






Vi G S" 


o = E 

jeS 

(107) 

since the states i G S” couple to non-selected states only. 
In these equations the coefficients for the states j G 


S' are determined already by Fqs. (103). 


The statement that all subsets of selected states must 
form a single connected cluster (in the sense described 
above) can now be shown by noting that the assump¬ 
tion of several mutually unconnected clusters A, B, C, 
... leads to a contradiction. Let us denote the set of 
non-selected states directly coupled to the selected states 
of cluster X by 5^ and note that the mean particle 
current from one subset of non-selected states to 
another one S 2 is in leading order given by Js^Si — 

SjGSi + Rn^? - The to¬ 

tal current into S" then reads Jg,, = Jgng, = Js"S'^ + 
Jgiigi -)-•■■. It is directly given by summing the right- 


hand-sides of Fqs. (107). Consequently, it vanishes in the 


steady state as it should, Jgn = 0. The total current into 


cluster A reads Jg/ = Jg> g,, + J. 


^S\S‘ 


+ J. 


S\S' 


•. Ob¬ 
viously, it should also vanish in the steady state. How¬ 
ever, generically this is is not possible for more than a 
single cluster. Namely, (without fine tuning) the indi¬ 
vidual terms Jg'^gn, Jg'^g'^i ■•■j containing the coeffi¬ 
cients determined from Eqs. (103) and (107), can 


neither be expected to vanish individually nor to can¬ 
cel each other. In contrast, for a single cluster, one has 
^ 5 ^ = Js'j^S" — ~Js" — required. 

From the rather technical discussion of the preceding 
paragraphs, we can now draw several important conclu¬ 
sions. First of all, Eqs. (97) and (98) imply that Bose 


selection is still predicted to occur, i.e. only a subset S 
of the single-particle states have occupations that grow 
with the total particle number 


Tli = ViN. 


(108) 


Second, the asymptotic occupations of the non-selected 
states are still determined by the first-order coefficient 


so that their occupations saturate for large N. (In 
contrast, if r') would describe the leading contribution 
to the occupations of a state i, it would become unpop¬ 
ulated in the limit of large particle numbers). This is 
true also for states contained in S" that are not directly 
coupled to a selected state. Both conclusions, Bose selec¬ 
tion and saturation, are confirmed W the preasymptotic 
state that can be observed in Fig. |^for 10^ ^ ^ 
which is approximately given by the asymptotic state of 
the rate matrix [Eq. . 

Finally, a third conclusion is that for rate matrices that 
are not fully connected the set of selected states S is not 
determined by the conditions (80) anymore. Namely a 


negative fii guarantees a positive asymptotic occupation 
of a non-selected state i G S', but not for a non- 
selected state i G S". This implies that we cannot apply 
the efficient algorithm presented in subsection |IV G| in 
order to find the set of selected states (neither can the 
algorithm of reference m be used, which is also based 
on the conditions (80)). It seems likely that the set of 


selected states of the mean-field equations is still unique 
and determined by the requirement of having positive 
occupations, as the full many-body master equation pos¬ 
sesses a unique steady state. However, unlike in the case 
of fully connected rate matrices, we have no proof for this 
statement. 

Let us illustrate the above reasoning using the min¬ 
imal example given by the rate matrix defined in 
Eq. (95) of the previous subsection. The corresponding 


rate-asymmetry matrix reads 





(109) 


Thus, if we were allowed to solve the problem (80) to find 


the set of selected states and the asymptotic occupations, 
we would find two disconnected clusters of selected states 
given by = {1} and S 2 = {3}. Namely, 


= r 



0 

r{vi - 21 ^ 3 ) I (110) 

0 


solves problem (80) non-uniquely for 0 < < 2/3 

and z /3 = 1 — vi. However, this is not the true so¬ 


lution. Namely Eq. (105) for a = 1 simplifies to 0 = 

' 2 i(i? 2 i + ^ 21 * 22 ^^) = 1 ^ 1(2 + 1 ^ 2 ^'^) from which 1/3 = 0 fol¬ 
lows in contradiction to Eqs. (80). This demonstrates 


that Eq. (80) cannot be used in order to determine the 


selected states in the case of non fully connected rate 
matrices. 

From Fig. where i?“ describes the preasymptotic 
regime (10^ ^ ^ 10^)j one can infer that only state 3 

will be selected. Let us, therefore, solve Eqs. ([^, d^, 
( fTMt , and dio^ for the ansatz 

5 = {3}. 


(Ill) 
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FIG. 10. (color online) Example for Bose selection of two 
uncoupled states. Mean occupations obtained from mean- 
field theory (hi blue, h 2 green, ha red, hi orange) vs the 
total particle number N for the rate matrix (1151, which is 
visualized in the inset (line widths reflect rates). The two 
selected states 2 and 4 are not coupled directly. 


The zeroth order equations (971 and (981 are solved triv¬ 
ially by 


Vs = 1 , 


Vi = 1/2 = 0 . 


( 112 ) 


Then ^2^^ is obtained from Eq. ([103| and reads 


2 Aa 

^23 


TDa 

^23 ^ ^ 


(113) 


while Eq. (1051 is trivially fu lfilled since is empty. 
Finally, results from Eq. (106) for 1 = 1, 


,(i) 


Da 

^ 12^2 


rya _ /la ,, 

-^21 ^ 12^2 


(1) 


1 

3’ 


(114) 


We can see that the initial assumption S = {3} is con¬ 
firmed by the fact that we obtained meaningful positive 
occupation numbers. The just-obtained asymptotic oc¬ 
cupations for the rate matrix are plotted as dotted 
lines in Fig. and provide a good description of the 
preasymptotic state. 

For completeness, we will finally present a simple ex¬ 
ample for a situation where the set of selected states con¬ 
sists of two uncoupled subsets. It is given by a model of 
four states with rate matrix 


/ 0 2 0 1 \ 
10 3 0 
0 10 4 
V 5 0 1 0 / 


(115) 


The occupations plotted in Fig. show that the set of 
selected states contains the two uncoupled states 2 and 
4, 


iS = 5iU52 , with = {2} and ^2 = {4}. (116) 


FIG. 11. (color online) Comparison between the dynamics 
of the mean-field equations (311 [solid lines] and the Lotka- 
Volterra equations (1171 [dashed lines] for the rate matrix 
(951, which is visualized in the inset (line widths reflect rates), 
for N — 300 particles initially uniformly distributed. While 
in both cases the population in state 2 (green) decays on an 
intermediate time scale, the population in state 1 (blue) de¬ 
cays on a longer time scale in the mean-field equations only 
leading to a single condensate in state 3 (red). 


The case of zero rates has recently also been discussed 
by Knebel et al. for the Lotka-Volterra equations of 
motion [151 HT] . 


ni = ni'^ AijUj. 


(117) 


These equations correspond to the leading-order high- 
density approximation (591 of the mean-field equation 


(31), with (T = 1 for bosons. These leading-order equa¬ 


tions describe the dynamics of the Bose gas on an in¬ 
termediate time scale, before, eventually the sub-leading 
terms of Eq. (31), which are linear in the occupations, be¬ 


come relevant and determine the steady state. Knebel et 
al. show that under the evolution described by Eq. (117) 


the occupations of some states i die out exponentially 
fast, while the other states retain non-zero occupations. 
Interestingly, those states retaining non-zero occupations 
are determined by the very same condition (80) that we 


found to determine the selected states for fully connected 
rate matrices. That means in the case of fully connected 
rate matrices the selected states are determined already 
by the leading-order equation (117). Note the general dif¬ 


ference between the mean-field equation on the one hand 
and the Lotka-Volterra equation on the other. While in 
the first case the non-selected states retain a small but 
non-zero occupation, they die out completely in the latter 
case. 

In the case of non fully connected rate matrices 
Eqs. (80) still determine uniquely which occupations die 


out under the dynamics of Eq. (117) m- Thus, the con¬ 
ditions (80) still describe a dynamical selection mecha¬ 


nism happening on an intermediate time scale. How¬ 
ever, in order to compute the (true) steady state ap¬ 
proached in the long-time limit, also higher-order equa- 
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tions [Eqs. (103), (104) and (106)] have to be taken into 
account. As a result, the set of selected states in the 
steady state can be different from that obtained from 
conditions (80). 

Let us illustrate the above reasoning using the exam¬ 


ple of the not fully connected rate matrix (95). Fig. 11 


shows the different dynamics of this system for both the 
full mean-field equations © and for the Lotka-Volterra 
equations of motion (59). In the limit of large N the pop¬ 
ulation of state 2 decays on the intermediate time scale, 
because the conditions (80) predict a extinction of occu¬ 
pations 712 on the level of the Lotka-Volterra equations 
[see Eq. (110)]. Eventually, however, when higher-order 
terms become relevant in the full mean-field equations 
of motion, also the occupation of state 1 d ecay s so that 
only state 3 is selected as predicted in Eq. (Ill). This is 
contrasted by the Lotka-Volterra system, which remains 
in the situation with two condensates in the states 1 and 
3. 


J. Asymptotic theory beyond mean field 


Our theoretical description of Bose selection has so far 
been based on mean-field theory. The data presented in 
Fig.i for a tight-binding chain in and out of equilibrium 
suggests that mean-field theory provides a rather good 
approximation to the mean occupations. Namely, devi¬ 
ations between the mean-field results (thick solid lines) 
and the exact Monte Carlo data (crosses) are visible only 
for non-selected states. And where visible deviations oc¬ 
cur they are still rather small and captured by the aug¬ 
mented mean-field theory (thin solid lines) introduced in 


Section 111 C Such good agreement can generally not be 
expected for the number fluctuations of macroscopically 
occupied selected modes, since mean-field theory does 
not comply with the conservation of the total particle 
number. 

In this subsection we will investigate corrections to 
mean-field theory in the asymptotic limit of large total 
particle number N, as they are described by the aug¬ 
mented mean-field theory. For simplicity, we will con¬ 
sider the case of fully connected rate matrices. We will 
explain why mean-field theory accurately describes the 
occupations of the selected states and that their correla¬ 
tions, such as number fluctuations, deviate from mean- 
field theory in a universal fashion. Moreover, we will 
argue that the set of selected states is well described by 
a Gaussian state projected to the space of sharp particle 
number N. 

Within the augmented mean-field theory (Section 


IIIC) the state of the system is described not only by 


the mean occupations = {fii), like in mean-field the¬ 
ory, but also in terms of the non-trivial two-particle cor¬ 
relations Qj = {fiifij) — In order to derive an 

augmented mean-field theory for the asymptotic limit of 
large total particle numbers N, we do not only expand 
the mean occupations with respect to the inverse par¬ 


ticle number, but at the same time also the non-trivial 
two-particle correlations, 

n, = Nv, + + • • • , (118) 

a. = • • • • (119) 

Moreover, we choose again the normalization conditions 
^zz, = l, (120) 

i i 

which fix the total particle number N = fii in leading 
order, as well as the conditions 

= ( 121 ) 

ij ij 


ensuring that the fluctuations of the total particle num- 

ber AN = Cij vanish. _ _ 

We now insert the expansions (118) and (119) into 
the augmented mean-field equations (45) and (46), with 
cr = 1 for bosons and with the left-hand side set to zero 
in order to obtain the steady state. In the resulting equa¬ 
tions we ask that all terms belonging to a certain power 
of N vanish independently. In this way, we obtain the 
set of coupled non-linear equations 


0 — 'y ] Aij [I'iVj + ^ij ], (122) 

3 

0 — y ( -\- AijUi^j^j (A/^j F Aij^i3j^i^i\ . 

(123) 


for the leading order. 

Remarkably, we can solve these equations by making 
the simple ansatz 


^ki — k^kk'i) ( 124 ) 


for the leading non-trivial correlations , with x being a 
free parameter. The relative weight of both terms in the 
bracket is chosen such that the condition (121) is obeyed. 
Entering the ansatz (124) into Eqs. (122) and (123) re¬ 


duces these equations to the much simpler conditions 


V, E = 0- 


(125) 


These equations are identical to the leading-order con¬ 
ditions ( |7l| ) of the asymptotic mean-field theory. Using 
the same arguments as in the conventional asymptotic 
mean-field theory, we have to conclude that the solution 
must be of the form 


E! “ 9, 

i € (S, 

(126) 

jes 



Vi = 0, 

iiS. 

(127) 

































25 


Equations (126) and (127) imply Bose selection. Only 
a subset S of selected states have non-vanishing occu¬ 
pations in leading order. The set S has to be deter¬ 
mined by the requirement to have positive occupations 
both for selected and non-selected states. The asymp¬ 
totic occupations of the latter are given by and have 
to be determined in the next order. Note that the set S 
obtained within the augmented theory can be different 
from the one obtained within mean-field theory. Namely, 
the occupations of the non-selected states differ in both 
theories, so that in the augmented theory, e.g., a tran¬ 
sition where S changes might be shifted away from the 
critical mean-field parameter. However, as long as the 
set of selected states is the same in both theories, the 
mean-field result for the asymptotic occupations of the 
selected states is not corrected anymore. This explains 
the excellent agreement between mean-field theory, aug¬ 
mented mean-field theory, and Monte-Carlo results for 
the selected-state occupati ons in Fig.[^ 

According to the ansatz (124), we find the asymptotic 
correlations among the selected states to be given by 


{fiifij) = (1 — x)ninj + xn^dij, i,j G S. (128) 

This is an intriguing result. It implies that the cor¬ 
relations and fluctuations are determined solely by the 
mean occupations and a single parameter x. The scaled 
two-particle correlations for particles in different selected 
states, 

= 1 - a;, (129) 

asymptotically approach all the same value, which is re¬ 
duced by X with respect to the mean-field result. 

This very same parameter x also determines the 
asymptotic number fluctuations of the Bose selected 
modes, 

Anf = Cm = xN‘^{ 1 — Vi)vi = x{N — ni)fii, i G S. (130) 


This equation implies that (in leading order) the number 
fluctuations vanish if we have a single condensate in the 
state i = k, so that = 1. This is a consequence of the 
conservation of the total particle number that is incorpo¬ 
rated in the augmented mean-field theory. It contrasts 
with the Gaussian result (35) obtained within the non- 
number-conserving mean-field theory, which for bosons 
(tj = I) reads Anf = fii + nf = + \/N)vi. Note, 

however, that as soon as a system features several con¬ 
densates (macroscopi cally occupied selected states), their 
number fluctuations (130) will typically be of the order 
of the total particle number. This reflects the fact that 
each condensate is effectively in contact with a particle 
reservoir given by the other ones. 

The requirement Anf > 0 tells us that x is positive. 
Moreover, it is reasonable to assume that the number 
fluctuations will not be much larger than those obtained 
within the non-number-conserving mean-field theory, so 


that Anf < h| for all i G S. Thus, an estimate for an 
upper bound for x is determined by the selected state i 
with the smallest occupation h, = ViN. Therefore, 


0 < a; < 


I - 


t'min = mint'* 
i£S 


(131) 


The precise value of x has to be obtained, however, from 
the first-order equations. These equations are rather in¬ 
volved and we will not discuss them here. They also de¬ 
scribe small beyond-mean-field corrections for the asymp¬ 
totic occupations, correlations, and fluctuations of the 
non-selected states. 


In Fig. 12 we compare the augmented theory (solid 


lines) with Monte-Carlo results (crosses with error bars), 
ordinary mean-field theory (dotted lines), and the asymp¬ 
totic prediction (128) for the selected states (dashed 
lines), using the model system of Fig. I^c). The com¬ 
parison with the Monte-Carlo data shows that the aug¬ 
mented mean-field theory provides an excellent approx¬ 
imation for the mean occupations Ui [panel (a)], where 
the ordinary mean-field theory shows small deviations 
for the non-selected states [see Fig. |^c)]. For the two- 
particle correlations {fiifij) shown in Fig. [T^b)-(f)], the 
augmented mean-field theory still provides a rather good 
description, though small systematic deviations with re¬ 
spect to the exact Monte-Carlo results are now visible, 
while mean-field theory is not reliable anymore. 

The relative number fluctuations An|/hf = 
for the selected states [panel (c)] show strong deviations 
from mean-field theory, once Bose selection sets in near 
N = 10^ [see panel (a)] so that the selected modes ac¬ 
quire “extensive” occupations. This agrees with our ex¬ 
pectation that mean-field theory is not able to describe 
the condensate fluctuations for a system with sharp par¬ 
ticle number. The condensate fluctuations are found to 


be consistent with the asymptotic prediction (130) for 
X « 0.018. Note that the selected state with the smallest 
occupation (roughly 4%) has asymptotic number fluctu¬ 
ations that are only half as large as the mean-field pre¬ 
diction, even though the other two condensates are large 
enough to serve as a reservoir. Thus x is roughly given 
by i^min/2 in agreement with the estimate (131). 

The other quantities displayed in Fig. are not ex¬ 
pected to exhibit such drastic deviations of orders of mag¬ 
nitude from mean-field theory, as we observed them for 
the condensate fluctuations. Panel (e) shows the scaled 
correlations (129) among the selected states. The aug¬ 


mented theory asymptotically approaches the universal 
value 1 — a;, with x « 0.018. Noticeable deviations of 
up to 30% occur before reaching the asymptotic regime, 
whereas the deviation from the mean-field result 1 be¬ 
come rather small asymptotically since a; <C 1. Similar 
behavior, i.e. larger deviations of up to a few tens of per¬ 
cent for small particle numbers that are reduced slightly 
in the asymptotic regime, can be observed also in the 
remaining plots of the figure. Panel (b) displays the rel¬ 
ative number fluctuations An| jnl for three exemplary 
non-selected states. Relative correlations pij between se¬ 
lected states and exemplary non-selected states as well 
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as among exemplary non-selected states are plotted in 
panel (d) and (f) respectively. 

A deeper understanding of the findings presented so far 
in this subsection, can be gained by noting that the se¬ 
lected states are asymptotically described by a projected 
Gaussian state. This can be seen as follows. For bosons 


in the steady state, the full many-body rate equation (13) 
takes the form 


o=E(i 


nj)ni 


RijPriji RjiPn 


(132) 


where Pn is the full occupation-number distribution. Let 
us accept that there is a group of selected states, whose 
occupations will grow with the total particle number N 
while all other occupations saturate. Asymptotically for 
Af —>■ 00 , we can then neglect all non-selected states and 
safely approximate {\ -\- rij) ^ rij, so that 


0 = rijUi 


RijPriji RjiPn 


(133) 


We can now show that this equation is solved by the 
projected Gaussian state (39). For this state one finds 
that pn-^ = e‘^'~^^pn for the probability of finding the 
system in the Fock state \nji) obtained from |n) by 
transferring one particle from i to j. Moreover, accord¬ 
ing to Eq. (33) one has e''* = 1 — l/fii ~ 1. Here we 
have used that, asymptotically, the mean occupations of 
the projected Gaussian state become identical to that of 
the non-projected Gaussian state. This implies the intu¬ 
itive statement that for the projected Gaussian state, the 
probabilities for finding the system in the almost identical 
Fock states and |n) asymptotically become identi¬ 
cal, Priji — Pn- Thus, plugging the proje cted Gaussian 
state into the right-hand side of Eq. (|133[), we obtain 


Pr, ’ 

i,jes 


R. .pVi-Vj -R 

J. Lij C- 


-Pr^Y. 


Tlj Tli ^ij 0, 


(134) 

since = —Aji. We have shown that asymptotically 
in the limit N ^ oo the full number distribution of the 
selected states is given by a projected Gaussian state. 
An important consequence is that mean-field theory pro¬ 
vides the exact asymptotic mean occupations of the se¬ 
lected states. Another consequence is that correlations 
{hifij) with i,j G S and, therefore, also the parameter a;, 
must be determined completely by the asymptotic mean 
occupations Nvi of the selected states. 


K. Heat flow through the system: the role of 
fragmented condensation and pseudotransitions 


Non-equilibrium steady states of a driven-dissipative 
quantum system typically feature a steady heat flow be¬ 
tween the system and its bath(s). This heat flow is de¬ 
scribed by Eq. 0 in the case of an autonomous system 


and by Eq. (10) for a periodically driven system. For 
bosons (cr = lyin a steady state, these equations read 

Qk = YiE^ - E,)Rf [{nf) + (n,n,)] (135) 


for the heat flow from an autonomous system into bath 
b and 

Q = - Ej - mhoj)R^E'> [(fi-) + {fiifij}] (136) 


for the heat flow from a Floquet system into a bath. In 
this subsection, we will investigate such heat flow in the 
regime of Bose selection. The dominant processes con¬ 
tributing to the heat flow will be identified. They are 
found to be given by transitions between different se¬ 
lected states and, for the Floquet system, also by pseu¬ 
dotransitions [cor responding to terms with i = j and 
m ^ 0 in Eq.(136)] associated with a selected state. 

In Fig. [^we present data obtained for a tight-binding 
chain that is driven between two heat baths, one of posi¬ 
tive temperature and a population-inverted one modeled 
by a negative temperature. This system corresponds to 
the one of Fig. Qc), but with the particle number fixed 
and with the relative coupling between both baths, 72 / 71 , 
varied. In panel (a) we plot the mean occupations versus 
the parameter p = {1 + 71 / 72 )”^: which increases with 
72 / 71 . One can observe several transitions. For p = 0, 
where the system is only coupled to bath 5 = 1 , a sin¬ 
gle state (the ground state) is selected as indicated by a 
large occupation. This corresponds to equilibrium Bose 
condensation. At a critical coupling to the second bath, 
near p = 0.2, three states become selected. Increasing the 
coupling to the second bath further, various transitions 
occur, where the set of selected states changes. Even¬ 
tually, roughly from p = 0.75 on only the most excited 
state will be selected, corresponding to the equilibrium 
situation at p = 1 , where the system is coupled to the 
population-inverted bath 2 only. In panel (b) we plot 
the heat flow from the hotter population-inverted bath 
through the system into the cooler positive temperature 
bath versus p. We can clearly see that the heat flow 
increases dramatically (by more than two orders of mag¬ 
nitude), when fragmented Bose condensation with more 
than just one selected state occurs. 

This effect, which has been reported already in ref¬ 
erence [23j . can be understood intuitively. Namely, in 
order to exchange energy with the system, the bath has 
to drive transitions between states i and j in the sys¬ 
tem. The larger the occupations of i and j, the larger 
will be the rate of the corresponding transition. There¬ 
fore, the most effective way of exchanging energy with 
the system is to drive transitions between two largely 
occupied states. And this is possible only if more than 
just one state is selected. This effect might be employed 
to control the heat conductivity of a bosonic system by 
switching between one and three selected states. 

In Fig. we show results for a periodically driven 
tight-binding chain coupled to a heat bath. This system 













27 




10^ 10^ 10^ lO'* ^ 10^ 


FIG. 12. (color online) Augmented mean-field theory (solid lines) versus mean-field theory (dotted lines) and Monte-Carlo 
simulations (crosses with error bars) for the tight-binding chain with parameters as in Fig. |^c). All colors are consistent 
with panel (a) [and also Fig. lines describing correlations between two states have alternating color, (a) Mean occupations 
[corresponding to thin solid ines in Fig. Be)]. (b) Relative number fluctuations An?/nf = (^njnl for three exemplary non- 
selected states, (c) Relative number fluctuations = (fulfil of the selected states, (d) Correlations gij between the 

selected and three exemplary non-selected states, (e) Correlations Qij among selected states i ^ j. (f) Correlations gij among 
exemplary non-selected states i ^ j- 


corresponds to the one of Fig. |^d), but with the particle 
number fixed and with the dimensionless driving strength 
varied. From the mean occupations plotted in panel 
(a), we can observe that for small 7 ^^ a single-particle 
Floquet state is selected, which is connected adiabatically 
to the ground state of the undriven system with 7 ^^ = 
0. Roughly at 7 ^^ = 0.25 and 7 ^^ = 1.5 the selected 
state changes in transitions, but still only a single state 
is selected. Only for a driving strength of about 7 (j = 2, a 
parameter window is reached, where three states become 
selected and acquire large occupations. 

The heat flow from the system into the bath is plotted 


in panel (b) of Fig. 14 In contrast to the autonomous 
chain, we can observe that the heat flow grows strongly, 
despite the fact that we have only a single selected state. 
This effect can be attributed to pseudotransitions [SS] 
associated with rates with i = j and m yf 0. In 

these processes the bath energy changes by mhoj, while 
the system’s state is not altered. Thus, the bath can 
effectively exchange energy with the system by driving 
pseudotransitions for a single strongly occupied (Bose se¬ 
lected) Floquet mode. This interpretation is supported 
by the dotted line, showing the share Q' of the heat flow 
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FIG. 13. (color online) Tight-binding chain with M — 20 
sites coupled to two heat baths. Parameters as in Fig. |^c), 
but for fixed N — 10^ and versus relative coupling strength 
(a) Mean occupations obtained from mean-field theory 
(solid lines), augmented mean-field theory (dashed line, indis¬ 
tinguishable from mean-field result) and Monte-Carlo simu¬ 
lations (crosses). Color code like in Fig. on the left-hand 
(right-hand) side the occupation decreases (increases) with 
energy, (b) Heat flow through the system from the hot¬ 
ter negative-temperature bath into the positive-temperature 
bath. 


FIG. 14. (color online) Periodically driven tight-binding chain 
with M = 20 sites coupled to a heat bath. Parameters as in 
Fig. I^d), but for fixed N — lO'* and versus dimensionless 
driving strength 7 ^^. (a) Mean occupations obtained from 

mean-held theory (solid lines) and Monte-Carlo simulations 
(crosses). Color code like in Fig. on the left-hand side 
the occupation decreases (increases) with energy, (b) Heat 
flow from the driven system into the bath obtained from 
mean-held theory (solid line), augmented mean-held theory 
(dashed line), and Monte-Carlo simulations (crosses). The 
dotted line is the mean-held heat how without the contribu¬ 
tion from pseudotransitions. 


not related to pseudotransitions, 


Q' = X! [(h,) -b {fiihj)]. 

(137) 

Away from the undriven limit = 0 and as long as only 
one Floquet mode i acquires a large occupation, Q' is 
typically two orders of magnitude smaller than the full 
heat flow and, thus, negligible. That means that prac¬ 
tically all the heat flow is based on pseudotransitions. 


the double sum in Eq. (1361 is dominated by the terms 


with i = j. Q' becomes significant only when several 
states have a large occupation. As one can clearly ob¬ 
serve in Fig. [T4|(b), this happens both near transitions, 
where two states are selected (see subsection IVF), and 
for 2 < 7 (j < 2.6, where three states are selected. Here an 
efficient heat exchange with the bath can be achieved by 
driving transitions between these largely occupied states, 
like for the autonomous system. 

In Fig. [T4|(b), we can also observe a noticeable differ¬ 
ence between the heat flow obtained from mean-held the¬ 
ory (solid line) and augmented mean-held theory (dashed 
line), in contrast to the autonomous system where both 
theories show very good agreement [on the logarithmic 
scale of Fig. mb) both lines overlap]. This is also a 
consequence of the strong impact of pseudotransitions 
in the condensate mode, which are determined by the 


condensate huctuations, a quantity that is overestimated 
by mean-held theory. This conhrms our conclusion that, 
thanks to pseudotransitions not present in autonomous 
systems, a bosonic Floquet system can be a good heat 
conductor even when most of its particles form a single 
Bose condensate. 

In conclusion, departing from equilibrium offers inter¬ 
esting possibilities to control the heat conductivity of a 
bosonic quantum system, which might be relevant for 
technological applications. 


V. IDEAL FERMI GASES 

In this section, we will briehy demonstrate that the 
theory of section jTT] and the methods presented in sec¬ 
tion m can also be employed to describe the properties 
of ideal Fermi gases. As a motivation, we note that the 
physics of such driven-dissipative Fermi gases will have 
to play an important role, for example, for the realiza¬ 
tion of Floquet topological insulators. These systems 
are based on lattice potentials that are forced periodi¬ 
cally in time such that they possess a topologically non¬ 
trivial quasienergy band structure giving rise to a quan¬ 
tized (spin) Hall conductivity, when one band is filled 
completely. Proposals for Floquet topological insulators 
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FIG. 15. (color online) Mean occupations versus total number 
of fermions N for a driven tight-binding chain with tunneling 
parameter J and M = 10 sites. The chain is coupled to 
a heat bath of temperature T = J at the first site and it 
is driven away from equilibrium by a time-periodic potential 
modulation at the last site of frequency huj = 1.5 J and driving 
strength yui = 2.3. Data obtained from mean-field theory 
(solid lines), augmented mean-field theory (dashed lines), and 
exact solution of the many-body rate equation (crosses). 


consider irradiated electronic systems like graphene [7] 
and semiconductor heterostructures [S3]; conceptually 
different schemes for the Floquet engineering of topo¬ 
logical band structures have been, moreover, proposed 
in the context of ultracold atomic quantum gases in 
optical lattices [m |72|. First experimental evidence 
of a (quantized) Hall conductivity in such systems has 
been observed with ultracold atoms in optical lattices 
mm- These systems are well isolated from their en¬ 
vironment. However, achieving this goal in an electronic 
solid-state systems, which cannot be viewed as isolated, 
is rather challenging. Namely, it cannot be expected 
that the periodically driven system in contact with the 
heat bath (given among others by phonons) will sim¬ 
ply form a band-insulating state with one band filled 
completely. Thus, one either has to resort to bath en¬ 
gineering in order to enforce a band insulating state 
mm or explore novel opportunities of tailoring inter¬ 
esting system properties related to non-thermal occupa¬ 
tions of (quasi)energy bands. In this section we will not 
address the issue of Floquet topological insulators, but 
present simple examples that show how the general for¬ 
malism of sections [^ and HI can be applied to compute 
non-equilibrium steady states of driven-dissipative Fermi 
gases. 

In Fig. we plot the mean occupations of a periodi¬ 
cally driven tight-binding chain of M = 10 states that is 
coupled to a heat bath and occupied by N spinless (i.e. 
spin-polarized) non-interacting fermions. The state is 
trivial not only for zero filling {N/M = 0), but as a conse¬ 
quence of Pauli exclusion also for unit filling {N/M = 1), 
corresponding to zero filling of holes. For intermediate 
filling N/M we find occupation numbers whose exact val¬ 
ues [obtained from solving the many-body rate equation 
@] are well described by mean-field theory. Residual 


deviations of the mean-field theory are cured within the 
augmented mean-field theory (Section HIC). 

As another example, we have computed steady states 
of a fermionic tight-binding chain of M = 100 sites (see 
section HD) and half filling (N=M/2). In Fig. 16 we plot 


the mean occupations of the single-particle states i of the 
chain versus their energy Ei = —2Jcos(A:i), where ki is 
the wave number of state i. In panel (a) the equilibrium 
situation is shown, where the system is coupled to a sin¬ 
gle heat bath of intermediate temperature T = J. The 
non-equilibrium system coupled to two baths of different 
positive temperature Ti = J and T 2 = 0.5J shows quali¬ 
tatively similar behavior, as can be seen from panel (b). 
In both situations (a) and (b) the occupations decrease 
with increasing energy. In striking contrast, the occupa¬ 
tions depend in a non-monotonous fashion on the energy, 
when the second heat bath is population inverted and de¬ 
scribed by a negative temperature. This can be seen in 
panel (c) and (d). Moreover, the distribution of occupa¬ 
tions depends sensitively on the structure of the system- 
bath coupling. Depending on whether bath 1 is coupled 
to the first site [panel (c)] or to the third site [panel (d)] 
the occupation of the ground state assumes either a local 
minimum or a local maximum. Thus, like in the bosonic 
case, already the ideal Fermi gas offers many possibilities 
of dissipative state engineering far from equilibrium. Ex¬ 
ploring these possibilities is, however, beyond the scope 
of the present manuscript. 


VI. CONCLUSIONS AND OUTLOOK 


In this paper, we describe several aspects of non¬ 
equilibrium steady states of driven-dissipative ideal quan¬ 
tum gases. We focus on systems of sharp particle num¬ 
ber that are driven away from equilibrium either by the 
coupling to two heat baths of different temperature or 
by time-periodic driving in combination with the cou¬ 
pling to a heat bath. We describe analytical and numeri¬ 
cal methods for treating these systems within the frame¬ 
work of (Floquet-)Born-Markov theory and apply them 
both to bosonic and fermionic quantum gases. On that 
basis, we work out a theory of Bose selection, a non¬ 
equilibrium generalization of Bose condensation, where 
multiple states can acquire large occupations. Also the 
possibility of bath engineering in a fermionic lattice sys¬ 
tem is pointed out. Our results demonstrate that already 
ideal quantum gases give rise to intriguing and unex¬ 
pected behavior, when they are driven into a steady state 
far from equilibrium. In the future it will be interesting 
to find applications for dissipative quantum engineering, 
e.g., in order to control the heat conductivity of a system 
in a robust fashion. On a theoretical level, it will be in¬ 
teresting to extend the formalism to systems exchanging 
particles with their environment and to include the effect 
of interactions. 






















30 



FIG. 16. (color online) Mean-occupations of the single-particle energy eigenstates in a tight-binding chain of M = 100 sites 
occupied by A'’ = M/2 spinless (i.e. spin-polarized) fermions versus the energy (in units of the tunneling parameter J). Data 
obtained from mean-field theory (solid lines) and exact Monte-Carlo simulations (crosses), (a) Equilibrium situation where 
the chain is coupled to one bath of temperature T — IJ. (b) The chain is driven away from equilibrium by two heat baths of 
different positive temperature (Ti = J and T 2 = 0.5J), coupled to the first and the next to last site with 71 = 72 . (c) Same as 
in (b), but now the second bath is population inverted and described by the negative temperature T 2 = — J. (d) Like in (c), 
but now the first bath is coupled to the third site. 
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Appendix A: Many-body rate equation from 
Lindblad master equation 

Here we derive the equations of motion for the many- 
body occupation probabilities Pn = (n|p|n), based on 
the Markovian master equation with the Liouvillian, 
Eq. 0. Replacing the single-particle operators |i)(j| by 
their representation in Fock space ajaj the equations of 
motion for the diagonal elements of the density operator 


take the form, 

Pn{t) ={n\p{t)\n) 

M 

hi=i 

-^{n\{p{t),a]aialaj}\n)y (Al) 

For i = j both terms inside the bracket cancel each other. 
For i 7 ^ j, we have a^aijn) = ± nj)\nji) and 

ajaia|aj|n) = n_,(l±ni)|n), where the upper (lower) sign 
applies to bosons (fermions). Thus, the master equation 
simplifies to 

M 

Pni^') — ^ ^ ^ ^ ^i)Vn{^^')\ t 

i,j=^ 

M 

= ^ (1 ± nj)ni [R^jPn-^{t) - RjiPnit)] . (A2) 
hi=i 
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wherein riji = (ni,..., — 1,..., rij + 1,...) denotes the 

occupation numbers obtained from n by transferring one 
particle from i to j. We have not explicitly excluded the 
i = j terms, since they still cancel. The second line was 
obtained by exchanging i and j in the second term. 


Appendix B: Equations of motion for mean 
occupations 


The equations of motion for the mean occupations read 


dt 


nfc(t) =tr ( nfe—/5(t) 


■-y^RijtT(nkalajp{t)a\di 






(Bl) 


where we have employed Eq. Q with the jump operators 
given by Eq. (12). The first term of the sum can be 
written like 


tr =tr {nkd^jdid\djp{t) 

+ {Sik - Sjk)tT 




a'jQialajpit 


(B2) 


Here we have used the invariance of cyclic permutations 
under the trace as well as the relation 




— {^ik ^jk) • 


(B3) 


This relation is valid for particles of either statistics, as it 
can be obtained both by employing either the commuta¬ 
tion relation [ai,a|] = Sij for bosons or the anticommu¬ 
tation relation {di,a)j} = Sij for fermions. We can now 
use 


djdioldj = nj{l ± fii) =F Sijfn, (B4) 

with the upper (lower) sign referring to bosons 
(fermions), to arrive at 

^hfc(t) ='^RijiStk - 5jk)ti[nj{l ± ni)p{t)) 

M , 

^ ^ j ^kj (f) ^ (^fc^_;)(^)] 

i=i ^ 

- Rjk[nkit) ± {nknj)(t)]\. (B5) 


Appendix C: Mean occupation and correlation in 
projected Gaussian state 


state 

Pproj CC P]\[PgP]\[ ^ (d) 

with 

Pn= Y. (C 2 ) 

n\ Yi ni—N 

is a non-trivial problem. This is why, already in equilib¬ 
rium it is typically much easier to treat a system in the 
grand-canonical rather than in the canonical ensemble. 
In this appendix we describe a method for computing 
expectation values 

(i)^r=tr(i/iproj) (C3) 

n 

for projected Gaussian states numerically. Here the sum 
constrained to Fock states of total particle num¬ 
ber N and Z = denotes the parti¬ 

tion function. 

We will focus on the mean occupations 

1 ^ 

(^*)jv = (C4) 

n 

and the two-particle correlations 
1 ^ 

(hjhj)Ar = — ^Uinj-e (C5) 

n 

The first expectation value can be written as 

(C6) 

rii 

wherein 

= Y^"" (“ mni) (C7) 

{nktkfSfi 

is the partition function of fictitious system obtained by 
the original one by removing the states Sr and filling it 
with Nr particles only. 

The second expectation value reads 

N N-m 

{n^fi,)^=-Y Y (C8) 

m—O rij—O 

The remaining partition functions can be calculated by 
exploiting the recursion formula m 

1 ^ 

= (C9) 

k+l 


Calculating expectation values, like mean occupation 
or second order correlations, for the projected Gaussian 


This enables the numerical treatment of systems with 
several thousands particles on M = 10 states. 
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Appendix D: Equations of motion for two-particle 
correlations 


In this appendix we derive the equations of motion 


Here, the first term and the anticommutator in Eq. (Dl) 


for the two-particle correlations (hfchi)(t) [Eq. (421] and 
rewrite this as equations of motion for the non-trivial 


correlations Cfci = (CfcC*) = - nkUi [Eq. (|^[. To¬ 

gether with the equations of motion for the mean occu¬ 
pations ni{t), Eqs. (45), they build the set of equation for 
the augmented mean-field theory described in Sec. |III C| 
Hereby we close the hierarchy of equation by assuming 
the three-particle correlations to be trivial. For the sake 
of a simple notation we will suppress the time argument 
in the following. 

The exact equations of motion for {fikhi) are obtained 
from the many-body master equation in Lindblad form 
Eq.(0 by multiplying it by fikni from the left and taking 
the trace, 


dt 


{hkhi) =tr {ukhi^ 

= ^i?i,tr( 


nknid\djpa}-ai 


2 


|/5, ajdiajdj 


})■ 


Invoking cyclic permutation under the trace and using 
Eq. (B31 we regroup the operators as 


a^jdinkfii =nknid^jdi + 

\{^li -f {Slk djk^fii 

{Sli dji)(^6ik 


(D2) 


form a commutator, which vanishes under the trace, 
tr{p[nkfii,a}jdialdj]) = 0. Applying also the operator 
relation Eq. (B4) we arrive at 


dt 


(hfcnj = ^ Rijtj:{ [{Su - Sj,)nk -f {Sit - Sjk)ni 


3^ 


+ (<5h - Sji){6ik - 6jk)] [nj{l ± fii) =F 5jifij\p]. 

(D3) 


The term Sjihj vanishes in combination with each of the 
(5-prefactors, leaving 


= ^Rij {Su - Sji) {{nkfij) ± {fikfijhi)) 
3,1 

+ {Sik - Sjk) {{n^nj) ± {hiJijhi)) 


(D4) 


Evaluating one of the two sums we arrive at 


3 

+ E - RjkiniJik) + R^j (nkfij) - Rj^{nkni)) 

3 

d~Sik ^ ^ (Rkj (nj dz (njflk)) Rjk (nk dz (nknj))) 

3 

Rik {nk i (nkni)) Rki {ni dz (flink)) , 


which is identical to Eq. (42). We separate the num¬ 


ber operators fik into their mean part and their fluc¬ 
tuations ({k = fik ~ f>-k- With that, the correlations 


read (nkfii) = nkni d- Cki with the non-trivial correlation 

Cki — (CkCi) ^Cld (flkniflj) — (CkCiCj) d~ nkC,ij d~ niQkj d~ 

fijCki d-nknifij- Now Eq. (D5) can be rewritten as 
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i ^ (CiCfcCj) '^iCkj '^kCij ‘^jCik '^k’^i^j 

3 

-4” ^ ^ {^Rkj \^i^j “t“ Cij\ \P‘i^k “1“ O/c] {^k^j Ckj] ^ji [^k^i Cfc^) 

3 

^ki E H- Rkj^ [‘^j'^k Cjk\ i^kj^j 

3 

^ {Rki H“ Rik) ['^k‘^i H“ Cki] {Rik'^k H“ Rki'^i)- 


(D6) 


To obtain the equations of motion for the non-trivial cor¬ 
relations Cfci, we subtract 

^(nkfii) = nk'^ {Rij{nj{l ± n^) ± Qj) 


i ^ C^i)) 

+ ni ^ (i?fc^(nj (l ± fik) ± Cfej) 


di 




— Rjk{nk(l ± rij) ± ^fc^)) 
from Eq. ( |D6[ ), to obtain 

- ^ ^ d" -^ij){{CiCkCj) + ^jCik) 

3 

^kjklkCij d“ ^ijklii^kj 

d" ^ ^ RjkCki d“ ^ijCkj RjiCki] 

3 

+5ki ^ [ di (iijfc + iifej) {fljUk + Cjfc) 

3 

d“ {^kjklj Rjk'klk'jl 


(D7) 


d“ Rik} (j^kk^i d“ C/ci) 
{Rik^k d“ Rkikki)- 


(D8) 


Finally, neglecting non-trivial three-particle correlations, 
iCkCi Ci) ^ Oj Oil® arrives at the non-linear set of equa¬ 
tions (46), which defines together with Eqs. (451 the aug¬ 
mented mean-field theory. 


Appendix E: Parameter-dependent solution of the 
auxiliary matrix A{p) 

For the auxiliary rate-asymmetry matrix A(p) given by 
Eqs. (1^ and (1^ the problem (80) takes the form 


ikiip) =E +pBij)i>j{p) 


3 

with 


Pi > 0 and fli = 0 for i G S{p), 
Pi = 0 and fli < 0 for i ^ Sip)- 


(El) 


r 


Together with Eq. (88) restricting B to have a cross-like 
structure, this implies 

T +P^ikbj - pSkjbi)Pjip) =0, iGSip). 

3eS{p) 

(E2) 

Let us now show that, unless a transition occurs where 
the set of selected states Sip) changes, the solution Pip) 
varies, apart from a normalization factor, linearly with p 
as written in Eq. (91). 

For that purpose we decompose the solution Ptip) like 
hip) = pf^ + ^Piip), iG Sip), (E3) 
where P^^^ is defined to solve 
-,( 0 ) 


E = 0, 

iG5(p) 


i G Sip). 


(E4) 


These equations possess a solution, since Aij is a skew 
symmetric matrix acting in the odd-dimensional sub¬ 
space spanned by the selected states. However, the P-°^ 
can be negative, as Sip) contains the selected states for 
the matrix Aip) and not for A. 

We can now distinguish two cases. If the state k is 
not c onta ined in the set of sel ected states, k ^ 5(fc), 
Eqs. (E2) simply reduces to Eq. (E4), so that we find the 
trivial parameter dependence 


hip) = P, 


(0) 


i G Sip), 


(E5) 


which complies with Eq. (91). If the state k is contained 
in the set of selected states, k G 5(fc), it is convenient 
to discard the normalization condition X)iG5(p) hip) = 1 
for the moment, in favor of requiring 


Vkip) = P, 


(0) 
k ’ 


i.e. 


APkip) = 0 . 


(E6) 


(E7) 


Note that this requires also to fix the solution of the 
homogeneous equations (E4) such that P^^ > 0, which 
we can always do. With that, all other states in Sip) 
obey 

E A.jAPjip) = pp^^h^, iGSip)\{k}. (E8) 

3&S{p)\{k} 
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This set of inhomogeneous equations possesses a solution, 
since Aij is a skew-symmetric matrix acting in the even¬ 
dimensional subspace spanned by the states of S{p)\{k}, 
which has no eigenvalue zero without fine tuning. The 
solution Ai>j{p) will depend linearly on the parameter p. 
Therefore, one finds that the i'i{p) depend linearly on the 
parameter p. 


with normalization constant C{p) > 0. One finds 


C ^{p)= ^ \i'f^+c,p 

ies(p) 


1+P ^ Ci, 
iG5(p) 


(Ell) 


Ui{p) = vf '' + Cip, i G S{p). (E9) 


In order to restore the normalization condition 
SiG 5 (p) ^iip) = 1) can now re-define 


Mp) = C{p) 


-m 


CiP 


i e S{p). 


(ElO) 


where the second equality holds if we c hoose 
= 1, which we always can. Eq. (ElO) 


-( 0 ) 


E 

implies that Eq. (91) is fulfilled also \i k G S{k). 


ieS(p) i 
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